%\documentclass[unicode,bachelor]{seuthesis} % 本科
\documentclass[master]{seuthesis} % 硕士
% \documentclass[doctor]{seuthesis} % 博士
% \documentclass[engineering]{seuthesis} % 工程硕士
\usepackage{amsmath}
 % 这里是导言区

\begin{document}
\categorynumber{000} % 分类采用《中国图书资料分类法》
\UDC{000}            %《国际十进分类法UDC》的类号
\secretlevel{公开}    %学位论文密级分为"公开"、"内部"、"秘密"和"机密"四种
\studentid{090543}   %学号要完整，前面的零不能省略。
\title{分治混合算法在计算电磁学中的应用}{}{Applications of Divide and Conquer Algorithm in Computational Electromagnetics}{}
\author{李\quad{}友}{LI You}
\advisor{崔铁军}{教授}{CUI Tiejun}{Prof.}
\coadvisor{俞文明}{副教授}{YU Wenming}{Associate Prof.} % 没有
                                % 可以不填
% \degree{工学硕士} % 详细学位名称
\major[12em]{电磁场与微波技术}
\defenddate{答辩日期}
\authorizedate{学位授予日期}
\department{信息科学与工程学院}{State Key Laboratory of Millimeter Waves\\ Department of Information Science and Engineering}
\duration{2009.9—2012.1}
\address{河海院2楼}
\maketitle

\begin{abstract}{区域分解\quad  混合\quad  电大尺寸\quad  快速算法}
随着电磁场数值分析技术的发展，电大尺寸结构的辐射散射问题受到人们的广泛关注，现代电磁学的研究大多是围绕这个问题展开的。
对于飞机、舰船级别等大型平台上的天线结构分析，普通的精确数值方法受限于未知量的数目已无法发挥效用。
针对这类电大目标和电小目标共存的问题，本文研究了基于物理光学法近似的高低频混合算法，把电大尺寸载体上的电流效应耦合到结构较精细的区域，而后通过求解这些精细区域得到整个系统的响应。

对于超电大尺寸目标或尺寸相近的群场景问题而言，上述方法受限于计算机的内存限制，且往往不满足应用高频近似的前提条件。
本文将此类目标结构表面分成若干子区域并进行缓冲区扩展，用积分方程对扩展后的子区域求解，有效地抑制了子区域边界上的电流奇异性，只需要两三次的外迭代即可准确的求解原问题。
同时利用了各个子问题相对独立存在的性质，对算法进行了分布式并行化的处理，在保证正确性的前提下最大程度的提高了计算效率。

通过使用基于惠更斯等效面的等效源(EPA)法，可以将物体表面的未知量转移到其近场的空间中，在不降低计算精度的情况下大幅改善目标的矩阵条件数，有效的解决了精度和收敛性无法同时进行优化的矛盾。
本文的最后基于此种方法对高低频混合及区域分解法进行了加速，并提出了一种稳定有效的普适性快速算法框架，达到了简化软件设计流程的目的。
\end{abstract}

\begin{englishabstract}{domain decomposition;\quad hybrid;\quad electrically large;\quad fast algorithm}
With the development of numerical analysis in electromagnetic fields, the radiation and scattering problem of electrically large structures has attracted widespread attention, and most studies of modern electromagnetics are focused on this area.
For the analysis of antenna on aircraft, ships and other large platform-level structures, the general numerical method is limited by the number of unknowns, and has been unable to be effective.
For such cases of small and large targets co-existence, this thesis investigates a high-frequency hybrid algorithm based on physical optics approximation.
Couple the currents of electrically large object to the finer structures, and then obtain the overall system's response by solving these finer regions.

For over electrically large targets or groups of similar size targets, the above method is limited by computer memory limitations, and often does not meet the prerequisite of high-frequency approximation.
This paper divide the surface of the structure into several extended sub-regions, and then use integral equation to solve these regions, effectively inhibiting the currents' singularity on the boundary.
Two or three times of the outer iteration can accurately solve the original problem.
While taking advantage of a relatively independent of each sub-problems, distributed parallel processing is taken to achieve the maximum computational efficiency on condition that the correctness of the problem is guaranteed.

By using the equivalence principle algorithm based on Huygens' equivalent surface, the unknowns on the surface can be transferred to the near-field volume of the space.
Without compromissing accuracy, the condition number of the target matrix is improved, effectively solving the conflict which convergence and accuracy can not be optimized simultaneously.
At the end of this thesis, the hybrid method and the domain decomposition method are accelerated by EPA, and a stable and effective framework for univeral fast algorithm is outlined, achieving the purpose of simplifying the software design process.
\end{englishabstract}
%\listoffigures
%\listoftables
%\begin{terminology}
%  本论文专用术语的注释表
%\end{terminology}

\begin{Main} % 开始正文

\chapter{绪论(前言)}
\section{课题背景}
现代社会无处不存在着大量的电磁设备和系统，如电视机、微波炉、移动电话、互联网、卫星通讯系统、雷达系统、医疗成像系统等等。
毫无疑问，电磁现象已经对人类的生活产生了深远的影响。
随着对电磁场的应用越来越广泛，对电磁场的研究也越来越深入。
1864年，麦克斯韦(Maxwell)用一组数学方程(附录\ref{MaxwellEquation})概括了宏观电磁场的基本规律，而后很长的一段时间内，人们对电磁场问题的分析多采用解析方法或渐近方法\cite{THEF}，如Mie在1908年提出的Mie级数解是所有算法中精度最高的，往往被作为数值方法比较的标准。
然而事实上只有一些经典的几何形状或结构较为简单的问题才有可能得到严格的解析解，如无限大平面、半无限大平面、球体、椭球体、无限长圆柱等。
实际应用中常常碰到的复杂电磁结构，解析方法往往无能为力，就是半解析的方法(如模式匹配法)也只能在个别问题中略显身手。

二十世纪六十年代以后，电子计算机的出现和发展开创了采用数值方法研究电磁问题的新时代。
其中最早出现并获得广泛应用的是矩量法和Yee提出的时域有限差分法，而后不久，常用于力学领域的有限元法也被证明可应用于计算电磁学中。
与经典方法相比，数值方法在所能处理的问题范围和复杂程度上表现出巨大的优越性。
因而，数值方法被广泛地应用到电磁学领域的方方面面，如微波与毫米波通讯、雷达、精确制导、导航和地质勘探、军用飞机、坦克、舰艇的隐身和反隐身、雷达系统设计和雷达目标识别、车载天线、卫星天线以及微波毫米波集成电路的设计等。

经典的电磁场数值分析算法大致可分为两大类：第一类基于微分方程，如频域有限差分(FDFD)\cite{TDRMCBFDM}、时域有限差分(FDTD)\cite{Yee1966}、有限元法(FEM)\cite{ITFEA}等，该类方法具有较强的通用性，具有简单、直观、适用范围广的特点，特别适合复杂边界和复杂介质电磁问题的求解。
但由于微分方程本身的截断误差会随着计算区域的变大而累积，在处理三维开放域的特大电磁问题时，此类方法便会遇到瓶颈。
第二类基于积分方程，如边界元法(BEM)\cite{BJYYYXY}等，该类方法求解的是电流或磁流，而电流或磁流只出现在边界面上，通常来说目标的体积大于“边界”，故就未知量数目而言积分方程类方法要比微分方程法少了很多。
由于格林函数的引入，不需专门考虑辐射边界条件。
基于这些特点，频域矩量法被视为数值算法中的精确求解方法\cite{XDJSDCXJC}，本文中研究的所有算法也都基于频域矩量法。
然而，频域矩量法所产生的矩阵是满矩阵，迭代求解的存储量级是$O(N^2)$，计算量级是$O(N^2)$，从对内存的需求和计算量来看，矩量法对电大电磁问题的分析和计算是一个挑战。

为解决上述问题，众多学者提出了可以大幅度降低内存和求解时间的快速算法，如快速多极子方法(FMM)\cite{FMMSOTDIE}、多层快速多极子方法(MLFMA)\cite{MFMAFESBLCO}、基于快速傅立叶变换的自适应积分方法、FFT法\cite{STSEPUFFTWCGM}和预校正P-FFT法\cite{PFSOTVIEF3IDO}等。
在多层快速多极子算法中，计算复杂度仅为$O(N\log N)$，所需内存可降到$O(N\log N)$，这个两个复杂度几乎已经达到目前理论水平的极限。
然而，当未知量数目$N$非常大时，MLFMA对内存和CPU的需求依然过大，现有的计算机资源往往都难以满足需求。
\section{国内外研究现状}
随着计算机技术和电磁场数值计算的发展，电大尺寸问题受到广泛关注。对于电大尺寸的电磁散射和辐射问题，矩量法的求解时间和对计算机资源需求的快速增长成为其主要的应用瓶颈。
几何光学法(GO)，几何绕射理论(GTD)\cite{GTOD,Keller1985,Keller1957}，物理光学法(PO)，物理绕射理论(PTD)\cite{EDAP,Bouche1993}等高频方法具有计算快速、所需的计算机存储量少的优点，可广泛应用于分析各类复杂目标的电磁散射特性，但存在计算结果精度较低的问题。

在早期的研究中，基于物理概念的近似方法和全波分析法相结合产生的诸如GTD-MM、PO-MM、SBR-MM、GTD-FEM、PO-FDTD等高低频混合算法\cite{Murthy1986,Hodges1997,Kaye1985,Newman1988,Jakobus1995,Thiele1992}得到了广泛的应用，这些方法在保证计算精度的同时又提高了算法的效率。
但由于使用的高频近似方法要求电磁波的频率足够高，限制了此类方法的使用范围。

近年来随着区域分解法(DDM)\cite{Li2008,Li2009,Hu2010,Li2008a,Li2011,ANISOTTDEFIE,QYFJSF-PWFFCSZJXJS}的提出，许多国内外学者将其应用于解决电磁问题。DDM不直接求解原电大问题，而采用“分而治之”的策略，把原来很大的求解区域划分成若干个相对独立的子区域，将原问题的求解转化为在各子区域上分别进行求解，并利用某种数据交换条件在子区域之间传递信息，最后得到原电大区域上的解。
与高低频混合算法相比，此种方法对电波的频率没有要求，大大增强了它的普适性，如今已有基于该方法的成熟的商业软件面世。

由于边界交换条件的限制，传统的DDM通常只采用某种单一的精确算法，如MoM、FDTD、FEM等，其效率还有提升的空间。
Chew等人提出的基于惠更斯(Huygens)面的等效源法(EPA)\cite{Li2006,Lu1995,Chew1993,Li2007,WAFIIM}，在使用DDM的基础上，同时去掉了对单一算法的限制，能同时满足精度和效率的要求，成为最成熟的多区域方法之一。
\section{本文内容和论文安排}
本文主要的研究对象为单一精细结构、电大尺寸以及两者共存情况下的辐射或散射问题。

本文所需实现的算法为矩量法、PO-MM高低频混合算法、基于矩量法的DDM和EPA等效源法。
目标是将电大尺寸、复杂结构的电磁场问题分解为两个或者更多相同或相似的子问题，直到最后子问题可以简单的直接求解，原问题的解即子问题解的合并。
通常来说，一种算法的收敛性和精度总是一组相互对立的矛盾，本课题的创新点在于建立一种通用的快速算法的框架，它能同时保证收敛性和精度达到一个最优值。

本文为基础性研究，旨在通过编程实现上述算法，实践中体会各种算法的优劣性，为本课题组的后续研究起到抛砖引玉的作用。

针对上述目标，本文具体章节安排如下：

第一章为绪论，介绍了本文的课题背景，讨论了当前计算电磁学的发展，论述了精确算法、高频近似算法以及混合算法的发展与现状，引出本文的研究意义。

第二章详尽的阐述了经典的全波分析算法——矩量法（MoM）。从讲述矩量法的数学原理开始，充分论述了任意结构三维形状物体的矩量法求解形式，给出了具体的公式推导，最后讨论了该算法的计算复杂度和内存效率，为后文基于矩量法的混合作好铺垫。

第三章介绍了物理光学法（PO）在高频假设下的感应流表达形式，并研究了一种基于电流修正的高低频混合算法(PO-MM)的原理及实现，最后将其与单一矩量法的结果进行对比，详尽讨论了该算法的优劣和计算效率。

第四章介绍了区域分解法(DDM)的理念，首先从数学原理上证明了区域分解法的可行性及效率，而后引出了基于表面积分方程的重叠型区域分解法，提出了解决区域边界电流奇异性的手段及实现方式，最后通过众多的算例验证了其效率和精度。

第五章首先针对前两章中所述算法的不足提出了基于Huygens原理的等效源法(EPA)，通过严格的数学推导和数值算例证明了其正确性，而后基于该方法建立了一套快速算法的框架，在保证收敛性的同时不损失精度，并能应用于任意的电磁算法中，拓展了它的实用价值。

在结束语中，回顾了本文的工作，并就后续的研究工作以及预期目标谈了一些看法。
\chapter{矩量法}
简而言之，离散积分方程数学表达形式的离散化方法称为矩量法\cite{JSDCXYL}。
它由Harrington于1968年的著作《Field Computation by Moment Method》\cite{FCBMM}中提出，由于积分方程法自动满足辐射边界条件，矩量法的求解空间限制在导体面或等效面上，因而其尤为适合求解开域问题，如目标的散射和辐射问题。
本章将从矩量法的数学原理开始，逐步讨论如何使用此方法解决三维结构的电磁散射或辐射问题。
\section{矩量法的数学原理}
通常所需求解的非其次积分方程有如下的表达形式
\begin{equation}
\mathcal{L}(f)=g\label{mom1}
\end{equation}
其中$\mathcal{L}$为线性算子，$g$是已知的激励或源函数，$f$为待求的未知函数。
在$\mathcal{L}$ 的定义域中将$f$用一组基函数$\{f_1, f_2, f_3,\dots, f_n\}$展开，得到
\begin{equation}
f=\sum_{n=1}^{N} \alpha_n f_n\label{mom2}
\end{equation}
这里$\alpha_n$是待定的标量，$f_n$被称为展开函数或基函数。
将\ref{mom2}式代入\ref{mom1}式，再由算子$\mathcal{L}$的线性性质可得
\begin{equation}
\sum_{n=1}^N\alpha_n\mathcal{L}(f_n)=g
\end{equation}

为了求得系数$\alpha_n$，还要再定义一组线性无关的检验函数(试函数)$\{\omega_1,\omega_2,\omega_3,\dots,\omega_m\}$，并对每个$\omega_i$取式(2.3)的内积，则有
\begin{equation}
\sum_{n=1}^N\alpha_n\left<\omega_i,\mathcal{L}\left(f_n\right)\right>=\left<\omega_i,g\right>
\end{equation}
这里$i=1,2,3,\dots,m$。上述方程可以写成矩阵形式
\begin{equation}
\left[Z\right]\left[I\right]=\left[V\right]
\end{equation}
其中
\begin{equation}\label{momV}
\left[V\right]=\left[\left<\omega_1,g\right>,\left<\omega_2,g\right>,\left<\omega_3,g\right>,\dots, \left<\omega_m,g\right>\right]^T
\end{equation}
\begin{equation}
\left[I\right]=\left[\alpha_1,\alpha_2,\alpha_3,\dots,\alpha_n\right]^{T}
\end{equation}
\begin{equation}\label{momZ}
\left[Z\right]=\left[
\begin{array}{cccc}
\left<\omega_1,\mathcal{L}(f_1)\right> & \left<\omega_1,\mathcal{L}(f_2)\right> & \cdots & \left<\omega_1,\mathcal{L}(f_n)\right> \\
\left<\omega_2,\mathcal{L}(f_1)\right> & \left<\omega_2,\mathcal{L}(f_2)\right> & \cdots & \left<\omega_2,\mathcal{L}(f_n)\right>\\
\vdots & \vdots & \ddots & \vdots \\
\left<\omega_m,\mathcal{L}(f_1)\right> & \left<\omega_m,\mathcal{L}(f_2)\right> & \cdots & \left<\omega_m,\mathcal{L}(f_n)\right>  
\end{array}
\right]
\end{equation}

上述\ref{momV}式称为右边项，\ref{momZ}式称为阻抗矩阵，均为已知量，如果矩阵$\left[Z\right]$非奇异，那么其逆矩阵存在。这样$\left[I\right]$便可求出
\begin{equation}
\left[I\right]=\left[Z\right]^{-1}\left[V\right]\label{mom3}
\end{equation}
通过\ref{mom3}式可以得到未知函数$f$的解。上述便是矩量法的完整离散化过程。
\section{基函数与试函数}
在任何一个特定问题中，离散积分方程的第一步就是选取基函数$\{f_n\}$和试函数$\{\omega_m\}$。
其中$\{f_n\}$必须是线性无关且完备的，并且使得它们的叠加形式\ref{mom2}能够很好的逼近$f$。
基函数主要分为全域基函数和分域基函数。
顾名思义，全域基函数是定义在整个求解域上的函数，这类基函数有收敛快、计算量少的优点，例如在计算线天线时使用的脉冲基函数就属于全域基函数。
但由于对很多问题，尤其是二维和三维问题，这类基函数很难构造，因而得不到广泛地使用。
分域基函数是在$f$定义域的各个分域上才存在的，它的特点是简单灵活，可以广泛应用于二维和三维问题中，尤其是Rao等人于1982年提出的RWG基函数\cite{Rao1982}几乎可以运用于任意形状的表面，本文的大部分算法也都基于这种基函数(见附录\ref{RWG})。

试函数的选择决定了测试的方式，最常见的有点匹配、线匹配和Galerkin匹配。
点匹配使用取样函数$\delta$做检验函数，在做内积后得到一系列特定离散点上的方程式；线匹配法需选取定义在连接两个长方形(或三角形)中心直线上的线性函数；Galerkin匹配是用基函数本身做检验函数的方法。
在计算难度上，点匹配最简单，线匹配次之，Galerkin难度最大，然而计算效果上以Galerkin匹配最为稳定，精度也最高，本文也将采用这种方法。
\section{表面积分方程}
自由空间中有一边界为$S$的理想导电体(PEC)，如图\ref{figure:PEC}所示。
\begin{figure}[htbp]
\centering
\includegraphics{PEC.eps}
\caption{PEC在平面波照射下}
\label{figure:PEC}
\end{figure}
在入射波$\bar{E}^{inc}$、$\bar{H}^{inc}$的照射下产生感应电流$\bar{J}$。由等效原理可知，散射场可等效为$S$上的感应电流$\bar{J}$在自由空间中产生的场，
\begin{equation}
\bar{E}^{sca}\left(\bar{r}\right)=
-j\omega \bar{A}\left(\bar{r}\right)-\nabla\Phi\left(\bar{r}\right)
\end{equation}
其中磁矢位$\bar{A}$定义为
\begin{equation}\label{A}
\bar{A}\left(\bar{r}\right)=\frac{\mu}{4\pi}\iint_{S'}dS' G\left(\bar{r},\bar{r}'\right)\bar{J}\left(\bar{r}'\right)
\end{equation}
电标位$\Phi$定义为
\begin{equation}\label{Phi}
\Phi\left(\bar{r}\right)=\frac{1}{4\pi\epsilon}\iint_{S'}dS' G\left(\bar{r},\bar{r}'\right)\sigma\left(\bar{r}'\right)
\end{equation}
$\epsilon$和$\mu$分别为自由空间中的介电常数和磁导率，$\bar{r}$和$\bar{r}'$分别表示场点和源点的位置。
面电荷密度$\sigma\left(\bar{r}\right)$通过电流连续性方程与面电流$\bar{J}\left(\bar{r}\right)$的散度相关联
\begin{equation}
\sigma\left(\bar{r}\right)=\frac{j}{\omega}\nabla\cdot \bar{J}\left(\bar{r}\right)
\end{equation}
$G\left(\bar{r},\bar{r}'\right)$为自由空间中的三维格林函数
\begin{equation}
G\left(\bar{r},\bar{r}'\right)=\frac{e^{-jkR}}{R},\quad R=\left|\bar{r}-\bar{r}'\right|
\end{equation}
因此式\ref{Phi}可改写为
\begin{equation}\label{newPhi}
\Phi\left(\bar{r}\right)=\frac{j}{4\pi\epsilon\omega}\iint_{S'}dS' G\left(\bar{r},\bar{r}'\right)\nabla '\cdot\bar{J}\left(\bar{r}'\right)
\end{equation}

由于导体表面满足电场切向分量为零的边界条件，即$\hat{n}\times\left(\bar{E}^{inc}+\bar{E}^{sca}\right)=0$，于是有
\begin{equation}\label{scattering}
\left[j\omega\bar{A}\left(\bar{r}\right)+\nabla\Phi\left(\bar{r}\right)\right]_{tan}=\left[\bar{E}^{inc}\left(\bar{r}\right)\right]_{tan}
\end{equation}
代入式\ref{A}和式\ref{newPhi}的结果可得到电场积分方程EFIE
\begin{equation}\label{EFIE}
\begin{split}
\bar{E}^{inc}\left(\bar{r}\right)&=\frac{j\omega\mu}{4\pi}\iint_{S'}dS'G\left(\bar{r},\bar{r}'\right)\bar{J}\left(\bar{r}'\right)\\
&+\frac{j}{4\pi\epsilon\omega}\nabla\iint_{S'}dS'G\left(\bar{r},\bar{r}'\right)\nabla '\cdot\bar{J}\left(\bar{r}'\right)
\end{split}
\end{equation}
同理由PEC表面的等效电流密度
\begin{equation}
\bar{J}\left(\bar{r}\right)=\hat{n}\times\bar{H}\left(\bar{r}\right),\quad \bar{r}\in S
\end{equation}
可得磁场积分方程MFIE
\begin{equation}\label{MFIE}
\bar{J}\left(\bar{r}\right)-\hat{n}\times\iint_{S}dS'\nabla\times G\left(\bar{r},\bar{r}'\right)\cdot\bar{J}\left(\bar{r}'\right)=\hat{n}\times\bar{H}^{inc}\left(\bar{r}\right), \quad \bar{r}\in S
\end{equation}
对EFIE和MFIE做线性组合可得到混合场积分方程CFIE
\begin{equation}\label{CFIE}
\textrm{CFIE}=\alpha \textrm{EFIE}+\left(1-\alpha\right)\eta \textrm{MFIE}
\end{equation}
式\ref{CFIE}中$\alpha$介于0～1之间，$\eta$为自由空间波阻抗$\sqrt{\mu/\epsilon}$。

事实表明，EFIE适用于分析任意结构的PEC散射体，而MFIE和CFIE只能用于封闭结构物体的分析。
MFIE产生的矩阵条件数最好但精度较低，EFIE产生的矩阵条件数差但精度较高，CFIE通过$\alpha$值的选择来兼顾收敛性和精度。
同时由于EFIE和MFIE的内谐振频率不同，因此采用CFIE分析封闭结构的散射体总是有效的。
\section{阻抗矩阵及右边激励项}
将物体的表面电流展开为
\begin{equation}
\bar{J}\left(\bar{r}\right)=\sum_{n=1}^{N}\alpha_n\bar{f}_n\left(\bar{r}\right)
\end{equation}
因此式\ref{A}和式\ref{newPhi}可改写成
\begin{equation}\label{discreteA}
\bar{A}\left(\bar{r}\right)=
\sum_{n=1}^{N}\alpha_n\frac{\mu}{4\pi}\iint_{S'}dS'G\left(\bar{r},\bar{r}'\right)\bar{f}_n\left(\bar{r}'\right)
=\sum_{n=1}^{N}\alpha_n\bar{A}_n\left(\bar{r}\right)
\end{equation}
\begin{equation}\label{discretePhi}
\Phi\left(\bar{r}\right)=
\sum_{n=1}^N\alpha_n\frac{j}{4\pi\epsilon\omega}\iint_{S'}dS'G\left(\bar{r},\bar{r}'\right)\nabla'\cdot\bar{f}_n\left(\bar{r}'\right)
=\sum_{n=1}^N\alpha_n\Phi_n\left(\bar{r}\right)
\end{equation}
其中
\begin{equation}
\bar{A}_n\left(\bar{r}\right)=
\frac{\mu}{4\pi}\iint_{S'}dS'G\left(\bar{r},\bar{r}'\right)\bar{f}_n\left(\bar{r}'\right)
\end{equation}
\begin{equation}
\Phi_n\left(\bar{r}\right)=
\frac{j}{4\pi\epsilon\omega}\iint_{S'}dS'G\left(\bar{r},\bar{r}'\right)\nabla'\cdot\bar{f}_n\left(\bar{r}'\right)
\end{equation}
用$\bar{f}_m\left(\bar{r}\right)$对式\ref{scattering}左右两边做测试可得
\begin{equation}
\iint_{S_m}\left[j\omega\bar{A}\left(\bar{r}\right)+\nabla\Phi\left(\bar{r}\right)\right]\cdot\bar{f}_m\left(\bar{r}\right)dS
=\iint_{S_m}\bar{E}^{inc}\left(\bar{r}\right)\cdot\bar{f}_m\left(\bar{r}\right)dS=V_m
\end{equation}
代入式\ref{discreteA}和式\ref{discretePhi}的结果
\begin{equation}
\sum_{n=1}^N\alpha_n\iint_{S_m}\left[j\omega\bar{A}_n\left(\bar{r}\right)+\nabla\Phi_n\left(\bar{r}\right)\right]\cdot\bar{f}_m\left(\bar{r}\right)dS=V_m
\end{equation}
以上是EFIE的矩量法解。输入阻抗矩阵为
\begin{equation}
Z_{mn}^E=\iint_{S_m}\left[j\omega\bar{A}_n\left(\bar{r}\right)+\nabla\Phi_n\left(\bar{r}\right)\right]\cdot\bar{f}_m\left(\bar{r}\right)dS
\end{equation}
右边激励项为
\begin{equation}
V_m^E=\iint_{S_m}\bar{E}^{inc}\left(\bar{r}\right)\cdot\bar{f}_m\left(\bar{r}\right)dS
\end{equation}

同理，MFIE矩量法解的输入阻抗矩阵为
\begin{equation}
Z_{mn}^H=\iint_{S_m}\left[\bar{f}_n\left(\bar{r}\right)-\frac{1}{2\pi}\bar{H}_n\left(\bar{r}\right)\right]\cdot\bar{f}_m\left(\bar{r}\right)dS
\end{equation}
其中
\begin{equation}
\bar{H}_n\left(\bar{r}\right)=\hat{n}\times\nabla\times\iint_{S_n}G\left(\bar{r},\bar{r}'\right)\bar{f}_n\left(\bar{r}'\right)dS'
\end{equation}
右边激励项为
\begin{equation}
V_m^H=2\iint_{S_m}\left[\hat{n}\times\bar{H}^{inc}\left(\bar{r}\right)\right]\cdot\bar{f}_m\left(\bar{r}\right)dS
\end{equation}
CFIE的输入阻抗矩阵和右边激励项均为EFIE、MFIE的线性组合：
\begin{equation}
Z_{mn}^C=\alpha Z_{mn}^E+\left(1-\alpha\right)\eta Z_{mn}^H
\end{equation}
\begin{equation}
V_m^C=\alpha V_m^E+\left(1-\alpha\right)\eta V_m^H
\end{equation}
\section{矩阵求解及远场计算}
通过本章前面的过程，我们将求解算子方程$\mathcal{L}\left(f\right)=g$转化为求解线性方程组，对线性方程组的求解有LU分解法、CG迭代法、GMRES迭代法等。
求解线性方程组之后，远区散射场可表示为
\begin{equation}
\bar{E}^{sca}\left(\bar{r}\right)=\frac{-j\omega\mu}{4\pi}\sum_{n=1}^N\alpha_n\iint_{S_n}G\left(\bar{r},\bar{r}'\right)\bar{f}_n\left(\bar{r}'\right)dS'
\end{equation}

对于远区，存在下列的近似
\begin{figure}[htbp]
\centering
\includegraphics[scale=0.7]{FFA.eps}
\caption{远场表达式中的近似}
\label{figure:FFA}
\end{figure}
\begin{equation}
\left|\bar{r}-\bar{r}'\right|\approx r-\hat{r}\cdot\bar{r}'
\end{equation}
因此散射场的表达式变为
\begin{equation}
\bar{E}^{sca}\left(\bar{r}\right)\approx\frac{-j\omega\mu}{4\pi}\sum_{n=1}^N\alpha_n\left[
\frac{l_n}{2A_n^{+}}\iint_{S_n^{+}}dS'\bar{\rho}_n^{+}e^{jk\hat{r}\cdot\bar{r}'}+
\frac{l_n}{2A_n^{-}}\iint_{S_n^{-}}dS'\bar{\rho}_n^{-}e^{jk\hat{r}\cdot\bar{r}'}
\right]\frac{e^{-jkr}}{r}
\end{equation}
对应的电场分量为
\begin{equation}
E_{\theta}^{sca}=\bar{E}^{sca}\cdot\hat{\theta}
\quad
E_{\phi}^{sca}=\bar{E}^{sca}\cdot\hat{\phi}
\end{equation}

为了更突出地反映物体的散射特征，通常使用均匀散射场能量来归一实际散射场能量，从而引入下面一个常常被使用的物理量：雷达散射界面(RCS)
\begin{equation}
RCS\left(\theta,\phi,\theta_i,\phi_i\right)=\lim_{r\to\infty}4\pi r^2\frac{\left|\bar{E}^{sca}\left(\theta,\phi\right)\right|^2}{\left|\bar{E}^{inc}\left(\theta_i,\phi_i\right)\right|^2}
\end{equation}
\section{数值算例}
\subsection{与Mie级数解比较}
算例1. 分别用三个积分方程来计算半径为一个波长的PEC球体(见图\ref{figure:PECsphere})的散射场，入射波频率为$f=0.3$GHz，波长为$\lambda=1$m。
球体表面用边长不大于$0.1\lambda$的三角面元剖分，有2048个单元，3072个RWG基函数，求解矩阵采用CG迭代法。
\begin{figure}[htbp]
\centering
\includegraphics[scale=0.5]{./data/sphere2.pdf}
\caption{半径一个波长的PEC金属球}
\label{figure:PECsphere}
\end{figure}

入射波方向为$\theta_i=0^\circ$，$\phi_i=0^\circ$，后图为计算$\theta=0$～$180^\circ$，$\phi=0^\circ$时的远场RCS与Mie级数解的对比，可以看出，采用三种方程都可以精确计算PEC散射体的RCS。
表\ref{table:IEcompare}列出了使用不同积分方程所消耗的时间，可以看出EFIE耗时最长，而MFIE和CFIE较短，验证了先前积分方程形式对矩阵条件数影响的推论。
\begin{table}[htbp]
\caption{三种积分方程耗时对比}
\label{table:IEcompare}
\centering
\begin{tabular}{ll}
\hline
积分方程 & 计算时间(s)\\
\hline
EFIE & 455.165\\
MFIE & 117.017\\
CFIE($\alpha$=0.2) & 89.443\\
\hline
\end{tabular}
\end{table}

\begin{figure}[htbp]
\centering
\includegraphics[scale=0.58]{./data/EFIE_VV.pdf}
\includegraphics[scale=0.58]{./data/EFIE_HH.pdf}
\newline
\caption{EFIE积分方程的结果}
\end{figure}

\begin{figure}[htbp]
\centering
\includegraphics[scale=0.58]{./data/MFIE_VV.pdf}
\includegraphics[scale=0.58]{./data/MFIE_HH.pdf}
\newline
\caption{MFIE积分方程的结果}
\end{figure}

\begin{figure}[htbp]
\centering
\includegraphics[scale=0.58]{./data/CFIE_VV.pdf}
\includegraphics[scale=0.58]{./data/CFIE_HH.pdf}
\newline
\caption{CFIE积分方程的结果}
\end{figure}
\subsection{与测量值比较}
算例2. 分别计算NASA标准体Almond在1.19GHz和Cone Sphere在869MHz时的单站RCS，并与实验测量值相比较。
由于这两种目标具有尖细结构，通过计算这种模型可验证矩量法在实际应用中的精度。
模型的摆放如图\ref{fig:almond_conesphere}所示，入射波均为$\theta_i=90^\circ$，$\phi_i=0$～$180^\circ$，计算结果如下。
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.45]{./data/almond.pdf}
  \includegraphics[scale=0.45]{./data/cone_sphere.pdf}
  \caption{NASA标准体Almond和Cone Sphere}
  \label{fig:almond_conesphere}
\end{figure}
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.57]{./data/Almond_1_19GHz_VV.pdf}
  \includegraphics[scale=0.57]{./data/Almond_1_19GHz_HH.pdf}
  \caption{Almond在1.19GHz时的单站RCS}
  \label{fig:almond_results}
\end{figure}
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.57]{./data/ConeSphere_869MHz_VV.pdf}
  \includegraphics[scale=0.57]{./data/ConeSphere_869MHz_HH.pdf}
  \caption{Cone Sphere在869MHz时的单站RCS}
  \label{fig:cone_sphere}
\end{figure}
\subsection{与商用仿真软件比较}
算例3. 由于本文后面大量使用商用仿真软件考量计算结果，因而此处将文中使用的矩量法和该仿真软件的计算精度进行对比。
同样计算一NASA标准体Ogive的散射场，Ogive的摆放如图\ref{fig:ogive}所示，入射波为$f=3.5$GHz，入射波的方向为$\theta_i=0^\circ$，$\phi_i=0^\circ$。
分别使用三种积分方程计算其$\theta=0$～$180^\circ$，$\phi=0^\circ$的远场RCS，并对比各种算法的余量误差和收敛速度。
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.43]{./data/ogive.pdf}
  \includegraphics[scale=0.58]{./data/Ogive_EFIE_VV.pdf}
  \caption{NASA标准体Ogive}
  \label{fig:ogive}
\end{figure}
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.57]{./data/Ogive_MFIE_VV.pdf}
  \includegraphics[scale=0.57]{./data/Ogive_CFIE_VV.pdf}
  \caption{NASA标准体及其散射场}
  \label{fig:ogive_RCS}
\end{figure}
\begin{table}[htbp]
\caption{三种积分方程计算Ogive的误差和收敛性}
\centering
\begin{tabular}{lll}
  \hline
  积分方程 & 迭代次数 & 余量误差\\
  \hline
  EFIE & 408 & 7.2185E-4\\
  MFIE & 57 & 0.0077\\
  CFIE & 68 & 0.0058\\
  \hline
\end{tabular}
\label{tab:Ogive}
\end{table}

上述结果再次验证了本文中使用的电场积分方程(EFIE)的精度高于磁场积分方程(MFIE)，但其阻抗矩阵条件数较差，收敛速度较慢；混合场积分方程(CFIE)结合了两者的优点，同时兼顾了计算精度和计算效率。
\section{本章小结}
本章从矩量法的数学原理出发，详尽阐述了RWG基函数的表现形式、表面积分方程的构建方法、矩阵元素的填充公式以及远场的计算方法，最后通过给出的数值算例证明了MoM的精度，并讨论了各种积分方程的优劣。本文的后续章节都将以本章的矩量法为基础，有针对性的研究电大问题的快速算法。
\chapter{物理光学法与矩量法的混合}
电磁散射的计算方法很多，大致可分为近似方法(或高频方法)和精确方法(或低频方法)两大类。
前一章节中介绍的矩量法属于低频方法，其具有原理简单易懂，适用于任何形状的物体且精确稳定的优点，然而待定的线性方程组系数矩阵是满阵，如用直接法求解，则计算机内存需$O(N^2)$，运算量达$O(N^2)$；如用迭代法求解，内存一样需$O(N^2)$，而每次迭代的运算量达$O(N^2)$(这里$N$为未知量的个数)\cite{JSDCXYL}。这样大的内存需求和运算量大大限制了矩量法的应用范围，使其仅仅适用于电小尺寸物体。

另一类的高频渐近方法具有计算效率高，对计算机内存和CPU速度要求较低的特点，但不足之处在于它对目标物体做了高频假设和无限大切平面近似，因此相对来说精度略低。
但当问题的电尺寸足够大时，高频方法所得到的结果往往已经能够满足人们对精度的要求，因此它通常被用于电大金属体电磁特性的预估。

针对上述数值方法性能各异、优缺点不一的特点，许多人提出了集各法之所长的解决方案，这便是电磁学中的混合算法。
本章将重点阐述物理光学法与矩量法混合的PO-MM方法。
\section{物理光学法}\label{section3.1}
物理光学法(PO)是由Macdonald于1912年提出的一种高频近似方法，通常用来计算电大尺寸的金属目标的散射问题。
在高频条件下，可以采用散射体表面的感应电流取代散射体本身作为散射场的源，然后对表面感应电流积分而求得散射场。

用PO法求感应电流时做了如下几点假设：
\begin{itemize}
\item 目标的特征尺寸远大于波长，即$D>>\lambda$；电磁波波数趋于无穷大，即$k\to\infty$
\item 散射体表面上只有被入射波直接照射的区域(照明区)才存在感应电流，不能被入射波直接照射的区域(阴影区)不存在感应电流，如图\ref{figure:PO}
\item 目标照明区上的感应电流特性和在入射点与表面相切的无穷大平面上的电流特性相同
\end{itemize}
\begin{figure}[htbp]
\centering
\includegraphics[scale=1.0]{PO.eps}
\caption{PO中的照明区与阴影区}
\label{figure:PO}
\end{figure}
假设表面$S$为一理想导体(PEC)面，根据等效原理有
\begin{equation}
\bar{M}^{sca}=-\hat{n}\times\bar{E}=0
\end{equation}
\begin{equation}
\bar{J}^{sca}=\hat{n}\times\bar{H}
\end{equation}
由于PO基于无限大切平面的假设，应用镜像原理可得
\begin{equation}
\bar{J}^{sca}=\hat{n}\times\bar{H}=2\hat{n}\times\bar{H}^{inc}
\end{equation}
因此PO区感应电流的强度可由入射磁场的强度直接求出，而不需要像矩量法般求解线性方程组。

下面通过一个算例验证PO方法的准确性，XY平面上放置一个边长为$5\lambda\times 5\lambda$的金属平板(图\ref{figure:5m_plate})，在入射波$f=0.3$GHz，$\theta_i=0^\circ$，$\phi_i=0^\circ$的照射下，对比单纯使用PO方法和快速多极子(MLFMA)法计算出的VV极化RCS。
\begin{figure}[htbp]
\centering
\includegraphics[scale=0.33]{./data/chapter3/5mPlate.jpg}
\includegraphics[scale=0.6]{./data/chapter3/POvsMLFMA.pdf}
\caption{$5\lambda\times 5\lambda$金属平板的散射}
\label{figure:5m_plate}
\end{figure}

可以看到，两种方法所得到的曲线趋势一致，在前向散射和后向散射角度附近吻合良好，而在$90^\circ$附近有一定误差。
这是由PO方法的近似性决定的，它假定表面连续且忽略阴影面以及边缘绕射等效应。
但对于电大尺寸的散射体，通常认为以牺牲一部分精度为代价而得到计算效率的极大提升是值得的。
\section{物理光学法与矩量法的混合}
鉴于PO方法有如上所述的特点，特别适合于矩量法混合，用于求解一个电小尺寸与电大尺寸共存的问题，例如导体面上的线天线、巨大金属平板上的球体等。
\subsection{区域的划分}
要使用PO-MM的混合方法，首先必须对整个计算目标进行区域划分，确定矩量法和物理光学法各自的应用区域。
划分的方案通常遵循下列规则：
\begin{itemize}
\item 精细复杂的结构通常归于MM区，光滑、平坦和电大结构通常属于PO区
\item 照明区一般为PO区，阴影区一般为MM区
\item 目标中的线状物体通常归于MM区，如图\ref{figure:diaplo_sphere}中的偶极子
\end{itemize}
上述准则仅作为以往经验的参考，可根据具体情况灵活更改。
\begin{figure}[htbp]
\centering
\includegraphics{diaplo_sphere.eps}
\caption{偶极天线置于导体球前方}
\label{figure:diaplo_sphere}
\end{figure}
\subsection{矩量法区的表达公式}
根据上一章阐述的矩量法可以得到线面混合结构的散射场表达式
\begin{equation}
\bar{E}^{sca}=\bar{\mathcal{L}}_J^E\bar{J}^{MM}+\bar{\mathcal{L}}_I^E I^{MM}
\end{equation}
\begin{equation}
\bar{H}^{sca}=\bar{\mathcal{L}}_J^H\bar{J}^{MM}+\bar{\mathcal{L}}_I^H I^{MM}
\end{equation}
其中$\bar{J}^{MM}$表示矩量法区的面电流，$I^{MM}$表示矩量法区的线电流。
$\bar{\mathcal{L}}_J^E$、$\bar{\mathcal{L}}_J^H$表示面电流产生电磁场的算子，$\bar{\mathcal{L}}_I^E$、$\bar{\mathcal{L}}_I^H$表示线电流产生电磁场的算子。

这四个算子的具体展开形式为
\begin{equation}
\begin{split}
\bar{\mathcal{L}}_J^E\bar{J}=&-\frac{j}{4\pi\epsilon\omega}\nabla\iint_{S'}\left(\nabla'\cdot\bar{J}\left(\bar{r}'\right)\right)\cdot G\left(\bar{r},\bar{r}'\right)dS'\\
&-j\omega\frac{\mu}{4\pi}\iint_{S'}\bar{J}\left(\bar{r}'\right)\cdot G\left(\bar{r},\bar{r}'\right)dS'
\end{split}
\end{equation}
\begin{equation}
\begin{split}
\bar{\mathcal{L}}_I^EI=&-\frac{j}{4\pi\epsilon\omega}\nabla\int_{L'}\frac{\partial I\left(\bar{r}'\right)}{\partial l'}\cdot G\left(\bar{r},\bar{r}'\right)dl'\\
&-j\omega\frac{\mu}{4\pi}\int_{L'}I\left(\bar{r}'\right)\cdot \hat{l'}\cdot G\left(\bar{r},\bar{r}'\right)dl'
\end{split}
\end{equation}
\begin{equation}
\bar{\mathcal{L}}_J^H\bar{J}=\frac{1}{4\pi}\nabla\times\iint_{S'}\bar{J}\left(\bar{r}'\right)\cdot G\left(\bar{r},\bar{r}'\right)dS'
\end{equation}
\begin{equation}
\bar{\mathcal{L}}_I^H I=\frac{1}{4\pi}\nabla\times\int_{L'}I\left(\bar{r}'\right)\cdot \hat{l'}\cdot G\left(\bar{r},\bar{r}'\right)dl'
\end{equation}
其中单位矢量$\hat{l'}$代表线上正电流的流向，$\nabla'$代表源点的面散度，$\nabla$代表场点的面散度。
矩量法区的未知电流同样用一组基函数$\bar{f}_n$、$g_n$和一组对应的复数系数$\alpha_n$、$\beta_n$来表示
\begin{equation}\label{SurfaceCurrent}
\bar{J}^{MM}=\sum_{n=1}^{N_J^{MM}}\alpha_n\cdot\bar{f}_n
\end{equation}
\begin{equation}\label{WireCurrent}
I^{MM}=\sum_{n=1}^{N_I^{MM}}\beta_n\cdot g_n
\end{equation}
$\bar{f}_n$和$g_n$通常分别使用RWG基函数和三角基函数，以保证面上和线上电流的连续性。
$N_J^{MM}$和$N_I^{MM}$分别代表面电流区和线电流区的未知量个数。

为了求解式\ref{SurfaceCurrent}和式\ref{WireCurrent}中的系数$\alpha_n$、$\beta_n$，可通过PEC表面的边界条件
\begin{equation}
\bar{E}_{tan}=0
\end{equation}
得到电场积分方程
\begin{equation}
\left(\bar{\mathcal{L}}_J^E\bar{J}^{MM}\right)_{tan}+\left(\bar{\mathcal{L}}_I^E I^{MM}\right)_{tan}=-\bar{E}_{tan}^{inc}
\end{equation}
代入式\ref{SurfaceCurrent}和式\ref{WireCurrent}的结果，并利用$\mathcal{L}$算子的线性特征
\begin{equation}
\sum_{n=1}^{N_J^{MM}}\alpha_n\cdot\left(\bar{\mathcal{L}}_J^E\bar{f}_n\right)_{tan}
+\sum_{n=1}^{N_I^{MM}}\beta_n\cdot\left(\bar{\mathcal{L}}_I^E g_n\right)_{tan}=-\bar{E}^{inc}_{tan}
\end{equation}
对上式使用Galerkin方法进行测试，可将积分方程转化为有$N_J^{MM}+N_I^{MM}$个未知量的线性方程组，由此可求得$\alpha_n$和$\beta_n$。
\subsection{物理光学法区的表达公式}
PO区的电流强度同样用RWG基函数进行展开
\begin{equation}\label{POcurrent}
\bar{J}^{PO}=\sum_{n=1}^{N_J^{PO}}\gamma_n\cdot\bar{f}_n
\end{equation}
与MM区不同的是，PO区的电流系数$\gamma_n$是通过\ref{section3.1}节中的等效原理直接求得的，而不是通过求解线性方程组。
式\ref{POcurrent}的拟合方式使得PO-MM混合法具有很强的灵活性，因为所有的目标表面都可以用三角面片进行剖分，而不必对PO区做特殊处理。
且RWG基函数的电流连续性也使得无论如何划分PO区和MM区，都能保证区域边界上的电流没有奇异性。
通过这样的PO近似，原线性系统中$N_J$的未知量可减少到$N_J-N_J^{PO}$，因此恰当的划分PO区和MM区可以有效的减小问题的规模。

根据\ref{section3.1}的结论，PO区的电流可以表示为
\begin{equation}\label{totalPOcurrent}
\begin{split}
\bar{J}^{PO}\left(\bar{r}\right)&=2\delta_i\cdot\hat{n}\times\bar{H}^{inc}\left(\bar{r}\right)\\
&+\sum_{n=1}^{N_J^{MM}}2\alpha_n\delta_{J,n}\cdot\hat{n}\times\bar{\mathcal{L}}_J^H\bar{f}_n\\
&+\sum_{n=1}^{N_I^{MM}}2\beta_n\delta_{I,n}\cdot\hat{n}\times\bar{\mathcal{L}}_I^H g_n
\end{split}
\end{equation}
式\ref{totalPOcurrent}右边第一项中$\bar{H}^{inc}$为入射波，$\hat{n}$代表物体表面的单位法向量，系数$\delta_i$代表阴影系数。
如果观察点$\bar{r}$位于入射波的阴影区，$\delta_i$为零，否则$\delta_i$为$\pm 1$，正负号取决于$\hat{n}$的定义。

传统的PO方法仅从入射波响应的角度求解表面感应电流，然而现在考虑这样一种情况：一个偶极天线(MM区)放置在一个反射面(PO区)前，入射波仅定义于偶极天线附近的一个小范围内，例如从中间点电压馈电。
此时，$\bar{H}^{inc}$就不存在于反射面上，意味着没有PO电流的产生，这显然是不正确的，因此矩量法区电流对PO区电流的影响必须被考虑进去。

式\ref{totalPOcurrent}右边项后两项中$\bar{f}_n$和$g_n$分别表示矩量法区面电流和线电流的基函数，$\alpha_n$和$\beta_n$分别表示各自的电流系数，$\bar{\mathcal{L}}_J^H$和$\bar{\mathcal{L}}_I^H$表示由电流产生磁场的算子，$2\hat{n}\times$得到PO电流，$\delta_{J,n}$和$\delta_{I,n}$依旧表示阴影系数。
因此这两项的含义就是，矩量法区电流产生的磁场在PO区产生的感应电流。

式\ref{POcurrent}中将PO电流用RWG基函数$\bar{f}_n$展开。
虽然$\bar{f}_n$不是一组正交基，但同样可以不通过求解线性方程组而得到系数$\gamma_n$的值。
这里引入两个新的单位矢量$\hat{t}_n^{\pm}$，它们位于第n条公共边的中点$\bar{r}_n=\frac{1}{2}\left(\bar{a}_{1,n}+\bar{a}_{2,n}\right)$，垂直于该条公共边，且分别位于三角面元$T_n^{\pm}$上，矢量方向如图\ref{figure:PO_RWG}。
\begin{figure}[htbp]
\centering
\includegraphics[scale=0.9]{PO_RWG.eps}
\caption{RWG基函数中各矢量示意图}
\label{figure:PO_RWG}
\end{figure}

基于$\bar{f}_n$垂直于公共边n的分量为单位矢量、而垂直于其它边的分量为零的特点，在第k条边的中点$\bar{r}_k$处可得到如下等式
\begin{equation}
\bar{f}_n\left(\bar{r}_k\right)\cdot\hat{t}_k^{\pm}=
\begin{cases}
1 & k=n\\
0 & k\ne n
\end{cases}
\end{equation}
因此对式\ref{POcurrent}左右两边同时乘以$\frac{1}{2}\left(\hat{t}_k^{+}+\hat{t}_k^{-}\right)$可得到
\begin{equation}
\gamma_k=\frac{1}{2}\left(\hat{t}_k^{+}+\hat{t}_k^{-}\right)\cdot\bar{J}^{PO}\left(\bar{r}_k\right)
\end{equation}
其中$k=1,2,\dots ,N_J^{PO}$。
将式\ref{totalPOcurrent}的结论代入上式中
\begin{equation}
\gamma_k=\tau_{i,k}+\sum_{n=1}^{N_J^{MM}}\alpha_n\cdot\tau_{J,n,k}+
\sum_{n=1}^{N_I^{MM}}\beta_n\cdot\tau_{I,n,k}
\end{equation}
其中
\begin{equation}
\begin{split}
\tau_{i,k}&=\left(\hat{t}_k^{+}+\hat{t}_k^{-}\right)\cdot\delta_i\cdot\hat{n}\times\bar{H}^{inc}\\
\tau_{J,n,k}&=\left(\hat{t}_k^++\hat{t}_k^-\right)\cdot\delta_{J,n}\cdot\hat{n}\times\bar{\mathcal{L}}_J^H\bar{f}_n\\
\tau_{I,n,k}&=\left(\hat{t}_k^++\hat{t}_k^-\right)\cdot\delta_{I,n}\cdot\hat{n}\times\bar{\mathcal{L}}_I^Hg_n
\end{split}
\end{equation}
\subsection{PO-MM混合算法表达式}
根据导体表面电场切向分量为零的边界条件可得到MM区的EFIE积分方程
\begin{equation}
\left(\bar{\mathcal{L}}_J^E\bar{J}^{MM}\right)_{tan}+
\left(\bar{\mathcal{L}}_I^EI^{MM}\right)_{tan}+
\left(\bar{\mathcal{L}}_J^E\bar{J}^{PO}\right)_{tan}=-\bar{E}_{tan}^{inc}
\end{equation}
注意与原MM区积分方程不同的是，此处考虑了2种不同的耦合效应：MM区电流对PO区的贡献，表现在$\bar{J}^{PO}$的表达式\ref{totalPOcurrent}中；PO区电流对MM区建立积分方程的贡献，表现为上式的左边第三项。
最后将各部分电流线性展开可得
\begin{equation}\label{POfinal}
\begin{split}
\sum_{n=1}^{N_J^{MM}}\alpha_n\cdot\left[\left(\bar{\mathcal{L}}_J^E\bar{f}_n\right)_{tan}+\sum_{k=1}^{N_J^{PO}}\tau_{J,n,k}\cdot\left(\bar{\mathcal{L}}_J^E\bar{f}_k\right)_{tan}\right]&\\
+\sum_{n=1}^{N_I^{MM}}\beta_n\cdot\left[\left(\bar{\mathcal{L}}_I^Eg_n\right)_{tan}+\sum_{k=1}^{N_J^{PO}}\tau_{I,n,k}\cdot\left(\bar{\mathcal{L}}_J^E\bar{f}_k\right)_{tan}\right]&\\
=-\bar{E}^{inc}_{tan}-\sum_{k=1}^{N_J^PO}\tau_{i,k}\cdot\left(\bar{\mathcal{L}}_J^E\bar{f}_k\right)_{tan}& \end{split}
\end{equation}
通过采用适当的检验函数，上式可转化为一个含有$N_J^{MM}+N_I^{MM}$个未知数的线性方程组，解这个线性方程组可求出矩量法区的电流系数$\alpha_n$和$\beta_n$，进而得到PO区的电流系数$\gamma_n$，整个系统的解也由此得到。

当计算目标中同时存在线电流和面电流时式\ref{POfinal}的表现形式还是比较复杂的，为了简化计算在求解线天线辐射问题时，通常将三维散射体全部归于PO区，MM区仅保留线天线。
由于MM区没有面电流的存在，且辐射问题仅在线天线上进行电压馈电，没有外加的激励波形，故式\ref{POfinal}中$\tau_{i,k}$，$\tau_{J,n,k}$，$\bar{\mathcal{L}}_J^H\bar{f}_n$均为零，公式可简化为
\begin{equation}
\sum_{n=1}^{N_I^{MM}}\beta_n\cdot\left[\left(\bar{\mathcal{L}}_I^Eg_n\right)_{tan}+\sum_{k=1}^{N_J^{PO}}\tau_{I,n,k}\cdot\left(\bar{\mathcal{L}}_J^E\bar{f}_k\right)_{tan}\right]
=-\bar{E}^{inc}_{tan}
\end{equation}
注意，在计算这类问题时涉及到线天线的矩量法，具体的公式形式见附录\ref{LineMoM}。

同理，在计算散射问题且目标中没有线状物体时，通常将精细结构归于MM区，剩余的电大尺寸部件归于PO区。
由于没有线电流的存在，式\ref{POfinal}中$\tau_{I,n,k}$，$\bar{\mathcal{L}}_I^E$为零，公式可简化为
\begin{equation}
\begin{split}
\sum_{n=1}^{N_J^{MM}}\alpha_n\cdot\left[\left(\bar{\mathcal{L}}_J^E\bar{f}_n\right)_{tan}+\sum_{k=1}^{N_J^{PO}}\tau_{J,n,k}\cdot\left(\bar{\mathcal{L}}_J^E\bar{f}_k\right)_{tan}\right]&\\
=-\bar{E}^{inc}_{tan}-\sum_{k=1}^{N_J^PO}\tau_{i,k}\cdot\left(\bar{\mathcal{L}}_J^E\bar{f}_k\right)_{tan}&
\end{split}
\end{equation}
\section{数值算例}
\subsection{MM区为线天线}
算例1. 考虑线天线与一个反射面放在一起的情况，下图为半波对称振子位于边长为$3\lambda\times 5\lambda$的矩形金属板前，天线轴垂直于平板(图\ref{figure:dipole_plate})，对称振子中间馈电，观察其$\phi=45^\circ$的方向图，并与快速多极子的结果对比(图\ref{fig:dipole_plate_result})。
\begin{figure}[htbp]
\centering
\includegraphics[scale=0.45]{./data/chapter3/dipole_plate.jpg}
\caption{偶极天线垂直于平板放置}
\label{figure:dipole_plate}
\end{figure}
\begin{figure}[htbp]
\centering
\includegraphics[scale=0.57]{./data/chapter3/dipole_plate_V.pdf}
\includegraphics[scale=0.57]{./data/chapter3/dipole_plate_H.pdf}
\caption{偶极天线位于平板附近的方向图}
\label{fig:dipole_plate_result}
\end{figure}

算例2. 现将PO区换成半径为$1\lambda$的金属球，对称振子放置方式如图\ref{figure:dipole_sphere}，观察$\theta=90^\circ$的垂直极化方向图(水平极化分量过小，不具参考价值)。
\begin{figure}[htbp]
\centering
\includegraphics[scale=0.5]{./data/chapter3/dipole_sphere.jpg}
\includegraphics[scale=0.58]{./data/chapter3/dipole_sphere.pdf}
\caption{偶极天线置于金属球附近}
\label{figure:dipole_sphere}
\end{figure}
\subsection{MM区为三维结构}
算例3. 半径为$1\lambda$的金属球置于$5\lambda\times 5\lambda$的平板上方，入射波为V极化，$\theta_i=0^\circ,\phi_i=0^\circ$(如图\ref{figure:sphere_plate})，观察其$\phi=0^\circ$的垂直极化RCS。
图中可以看出在$\theta=40$～$140^\circ$范围内精度有一定误差，而在前向和后向角度上吻合的却很好。
这也说明PO-MM混合有它算法自身的局限性，通常只用于计算目标的单站电磁特性。
\begin{figure}[htbp]
\centering
\includegraphics[scale=0.5]{./data/chapter3/sphere_plate.jpg}
\includegraphics[scale=0.6]{./data/chapter3/sphere_plate.pdf}
\caption{金属球置于平板上方}
\label{figure:sphere_plate}
\end{figure}

算例4. 考虑一个比较复杂的例子，f16战斗机的机身上有一金属天线，计算在频率为0.3GHz的入射波照射下的远场RCS。
f16战斗机的尺寸为16m$\times$10m，在这个问题中即为$16\lambda\times 10\lambda$的电尺寸。
为简化模型的描述，在这里可以不失一般性的用半径为0.5m的金属球体代表该天线。
机身除发动机和天线支架的部分为介质外，其余部分均为理想导体，模型中已将介质的部分抠去。
分别计算只使用PO和使用PO-MM混合情况下的结果，并与快速多极子(MLFMA)的结果进行对比(图\ref{fig:f16_sphere_results})。
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.45]{./data/chapter3/f16.pdf}
  \includegraphics[scale=0.45]{./data/chapter3/f16_sphere.pdf}
  \caption{f16战斗机的原始模型和带天线的简化模型}
  \label{fig:f16_Model}
\end{figure}
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.7]{./data/chapter3/f16_sphere_PO_MoM_MLFMA_VV.pdf}
  \caption{f16战斗机的雷达散射截面}
  \label{fig:f16_sphere_results}
\end{figure}

可以看到，单纯使用PO时，仅有前向和后向附近的曲线吻合的较好。
而使用PO-MM的混合后，$\theta=-100^\circ$附近的RCS也得到了修正，更接近精确值。
但和算例3一样，无论是哪种近似算法，在$\theta=40$～$140^\circ$范围内精度仍有较大误差。
\subsection{结果分析}
通过使用PO-MM混合方法，可以极大的减少目标未知量的数目。
上述例子中，半波长的线天线分成均匀的9段，面电流部分均按最长$0.1\lambda$的边长进行剖分，使用混合方法前后的未知量对比如表\ref{table:PO-MM_unknown}。
\begin{table}[htbp]
\caption{使用混合方法前后未知量数目的对比}
\label{table:PO-MM_unknown}
\centering
\begin{tabular}{llll}
\hline
算例 & 单一矩量法 & 使用PO-MM混合 & 未知数减少量\\
\hline
线天线+平板 &  5011 & 9 & 99.82\%\\
线天线+球 & 3011 & 9 & 99.70\%\\
球+平板 & 11288 & 3072 & 72.79\%\\
球+飞机 & 121843 & 3072 & 97.48\%\\
\hline
\end{tabular}
\end{table}

可以看出，使用了PO-MM混合算法后虽然精度上有所损失，但却在计算效率上得到了质的飞跃，因此特别适合被用于做电大尺寸模型的电磁特性估计。

我们再来考虑该方法在计算过程中的算法复杂度。式\ref{POfinal}中，$\tau_{J,n,k}$为$N_J^{PO}\times N_J^{MM}$的矩阵，算子$\bar{\mathcal{L}}_J^E\bar{f}_k$为$N_J^{MM}\times N_J^{PO}$的矩阵，$\tau_{J,n,k}\cdot\left(\bar{\mathcal{L}}_J^E\bar{f}_k\right)_{tan}$为两个矩阵的乘法运算，例如上述球+平板的例子中就需要求解一个$3072\times 8216$矩阵和$8216\times 3072$矩阵的乘积。
这部分的矩阵乘计算被认为是该算法中最耗时的部分，有的时候甚至使得混合方法的效率还不如普通矩量法，因此通常使用并行计算等技术对矩阵乘法进行优化和加速。后面的章节中还会介绍到另外一种降低矩阵规模的方法，通过这些手段，PO-MM混合方法的应用途径得到了极大的扩展。
\section{本章小结}
本章首先介绍了物理光学法的概念，紧接着在第二章矩量法的基础上，论述了一种高低频混合算法(PO-MM)在求解电大尺寸问题时的应用。
包括如何划分区域，如何离散PO区的电流，如何修正MM区的电流等关键性问题。
而后通过几个算例证明了该算法在求解辐射和散射问题时的有效性。
最后对比了使用混合算法前后的未知量数目，并分析了该算法的复杂度，提出了进一步提高算法效率的措施。
\chapter{区域分解矩量法}
上一章中所介绍的PO-MM混合方法可归类于一种电流修正的混合方法，它将高频区的未知电流耦合进矩量法区的积分方程中，而后通过求解一次矩量法区的电流系数即可得到整个系统的所有未知量。
然而正如上一章最后提到的，此类的混合方法需要进行大矩阵的乘法运算，某些程度上降低了算法的效率和适用范围。

本章将介绍另一类基于迭代的混合方法，此类方法同样将原问题分解为两个或者更多的区域，但每次仅计算单个区域的解，区域间的相互作用是通过多次迭代考虑进去的。
此类混合方法虽然需要对同一区域进行多次重复计算，但这些区域的规模通常较小，且省去了诸如PO-MM中的大矩阵相乘等耗时的步骤，不需要做高频假设和切平面近似，因此算法效率和易用性得到了进一步的提升。
\section{区域分解法的原理和数学描述}
在使用矩量法求解目标的未知电流时将积分方程转化为线性方程组，写成矩阵形式即
\begin{equation}
Ax=b
\end{equation}
其中激励矢量$b$被称为右边项(RHS)，$b$的长度取决于系统中未知量的数目。

求解上述矩阵方程有很多种方法，直接求解法是其中最简单易懂的，它通过矩阵$A$的逆得到
\begin{equation}
x=A^{-1}b
\end{equation}
然而这种方法的计算复杂度高达$O(N^3)$，当未知量的数目过大时便不再适用，取而代之的是迭代求解的方式。
迭代求解法将原矩阵方程修正为
\begin{equation}
MAx=Mb
\end{equation}
其中$M$为预处理矩阵，用于改善矩阵方程的条件数使其更快收敛，通常用CG、GMRES等迭代方法求解上述方程。

无论是传统的直接求解法还是迭代求解法，都只能在一台机器的内存上进行，而下面介绍的区域分解法将使得多机的分布式处理成为可能。

区域分解法(DDM)基于分治算法的原理，它不直接求解复杂庞大的问题，而是通过分区将其分解为两个或多个容易求解的子问题。
这种方法的关键之处在于如何通过边界条件的定义，使得分解后的子问题不破坏原目标中的电磁场连续性。
为了描述该方法的基本思想，此处将求解$Ax=b$的问题分解到两个子域
\begin{equation}\label{originalDDM}
\left[\begin{array}{cc}
A_{11} & A_{12}\\
A_{21} & A_{22}\\
\end{array}\right]
\left[\begin{array}{c}
x_1\\
x_2\\
\end{array}\right]
=\left[\begin{array}{c}
b_1\\
b_2\\
\end{array}\right]
\end{equation}
其中$A_{ii},x_i$和$b_i\left(i=1,2\right)$是区域$i$的系统矩阵、解向量和右边项，$A_{12},A_{21}$是两个区域之间的耦合矩阵。

为了并行求解各个子域的问题，使用一种常用的Jacobi型区域分解算法
\begin{equation}\label{Jacobi}
\left[\begin{array}{cc}
A_{11} & A_{12}\\
A_{21} & A_{22}\\
\end{array}\right]
=\left[\begin{array}{cc}
A_{11} & 0\\
0 & A_{22}\\
\end{array}\right]
+\left[\begin{array}{cc}
0 & A_{12}\\
A_{21} & 0\\
\end{array}\right]
\end{equation}
式\ref{originalDDM}使用迭代解法可变为
\begin{equation}\label{iterativeDDM}
\left[\begin{array}{cc}
A_{11} & 0\\
0 & A_{22}\\
\end{array}\right]
\left[\begin{array}{c}
x_1\\
x_2\\
\end{array}\right]^{(n)}
=\left[\begin{array}{c}
b_1\\
b_2\\
\end{array}\right]
-\left[\begin{array}{cc}
0 & A_{12}\\
A_{21} & 0\\
\end{array}\right]
\left[\begin{array}{c}
x_1\\
x_2\\
\end{array}\right]^{(n-1)}
\end{equation}
通过简单的代数运算，式\ref{iterativeDDM}转化为两个耦合的系统
\begin{equation}
\begin{split}
&{A_{11}x_1}^{(n)}=b_1-{A_{12}x_2}^{(n-1)}\\
&{A_{22}x_2}^{(n)}=b_2-{A_{21}x_1}^{(n-1)}
\end{split}
\end{equation}
设${x_1}^{(0)}=0,{x_2}^{(0)}=0$，可以得到初始值${x_i}^{(1)}={A_{ii}}^{-1}b_i$，它们分别为各个子区域的本地解。
这两个初始值通过耦合矩阵$A_{12}$和$A_{21}$的修正直到达到稳态。
\section{区域分解法的复杂度分析}
相对于原问题需要求解$A$矩阵来说，DDM仅需求解两个小得多的矩阵$A_{11}$和$A_{22}$，相比之下有着更高的内存使用效率。
假设原问题有$N$个未知数，通过DDM分解成$M$个同样大小的区域，如两个区域的情况下$N_1=N_2=N/2$，$N_1$和$N_2$为各子域的未知量数目。
原问题的内存需求量为$C_1N^\alpha$，其中$C_1$是常量而$\alpha>1$。
对于$M$个子域的分解方法，总的内存需求为
\begin{equation}
MC_1\left(\frac{N}{M}\right)^\alpha=M^{1-\alpha}C_1N^\alpha
\end{equation}
例如，当$\alpha=1.5$时，2个区域的DDM将节省30\%的总内存量。
另外还需要注意的是，使用DDM求解问题可以将内存进行分布式处理，进一步的增大了可求解问题的规模。

现在再来看时间复杂度，区域分解算法通常包括下列几个步骤：
\begin{enumerate}
\item 网格分区，此过程极快，几乎不占计算时间，且可预先处理好
\item 矩阵填充，矩阵$A$和$M$都可并行填充
\item 区域迭代，求解$MAx=Mb$
\end{enumerate}
直接求解系统方程$Ax=b$需要$C_2N^\beta$的时间复杂度，其中$C_2$为常量$\beta>1$(通常比$\alpha$大)。
而对于$M$个子域的DDM，主要的计算复杂度为
\begin{equation}
C_2\left(\frac{N}{M}\right)^\beta=m^{-\beta}C_2N^\beta
\end{equation}
因此，当$\beta>1$时，该方法呈现出超线性加速的特点。
\section{基于积分方程的区域分解法}
DDM方法有着极高的效率，因而首先被广泛应用于基于泛函的有限元(FEM)领域，然而近年来的研究证明它对基于积分方程的矩量法同样适用。

首先，任意形状的理想导体(PEC)表面被分成若干个子区域，如图\ref{figure:partition}实线包围的区域。
\begin{figure}[htbp]
\centering
\includegraphics[scale=0.6]{DDM_partition.eps}
\caption{矩量法区域分解示意图}
\label{figure:partition}
\end{figure}
若直接按上图中实线所示的分区求解，则得到的分块矩阵相互之间无公共元素。
从电流角度看，分区边界有强加的截断边界条件，不符合整个区域内部电流连续的分布规律。
这样子区域之间的迭代效率较低，通常会造成不收敛的情况。

为了保证收敛性和提高收敛速度，可采用重叠型的分区方法，即广义矩阵分解，模型上表现为子区域带有缓冲区，且与邻近区域有重叠。
为抑制边界电流的奇异性，Brennan等提出了前向后向缓冲迭代算法(FBBS)。
如图\ref{figure:subregion}所示，假设整个求解区域为$\Omega$，黑色区域表示第$i$个子区域$\Omega_i$，两片阴影区分别代表前向和后向的缓冲区$\Omega_{Fi},\Omega_{Bi}$。
\begin{figure}[htbp]
\centering
\includegraphics[scale=0.5]{subregion.jpg}
\caption{求解子区域与缓冲区}
\label{figure:subregion}
\end{figure}
当求解方向为向前迭代时，形成新的子区域$\Omega_i+\Omega_{Fi}$，同理向后迭代时形成$\Omega_i+\Omega_{Bi}$。
但每次解得扩展后的子域$i$的电流时，只保留$\Omega_i$上的电流，遗弃$\Omega_{Fi}$或$\Omega_{Bi}$上的电流。
被保留区域电流对其它区域的影响通过格林函数实现。

然而这种算法只是部分解决了子区域边界的电流奇异性，如前向迭代时只有前向边界上的奇异性被抑制，其它边界上扔存在奇异性，后向迭代亦然。
因此又有人提出对子区域的前向后向同时进行扩展，形成$\Omega_i+\Omega_{Fi}+\Omega_{Bi}$的新扩展域进行求解，同样仅保留$\Omega_i$上的电流。
这种方法也有其局限性，因为对于三维散射体，并不像无限薄平板一样可以定义其前向和后向。
但是却不难想到，既然三维散射体没有所谓的前后向存在，那索性将其所有边界均进行扩展。

因此原算法变为：每个子区域$\Omega_i$用其附近的缓冲区$\Omega_{bi}$进行扩展，依次对扩展后的区域$\Omega_i'=\Omega_i+\Omega_{bi}$求解(如图\ref{figure:partition}中虚线包围的区域)，舍弃缓冲区的电流，保留并更新原子区域$\Omega_i$的电流，这里定义各扩展子区域的补余区域$\bar{\Omega_i}=\Omega-\Omega_i'$。

根据上述原理，可建立子区域间的EFIE积分公式
\begin{equation}
\bar{K}_e\left(\bar{r},\bar{J}\right)=
jk\eta\left[\int_{\Omega_i}ds'\bar{G}\left(\bar{r},\bar{r}'\right)\cdot\bar{J}\left(\bar{r}'\right)\right]_{tan},\bar{r}'\in\bar{\Omega_i},\bar{r}\in\Omega_i'
\end{equation}
\begin{equation}
\bar{T}_e\left(\bar{r},\bar{J}\right)=
jk\eta\left[\int_{\Omega_i'}ds'\bar{G}\left(\bar{r},\bar{r}'\right)\cdot\bar{J}\left(\bar{r}'\right)\right]_{tan},\bar{r}'\in\Omega_i',\bar{r}\in\Omega_i'
\end{equation}
其中，$\bar{G}$为自由空间内的并矢格林函数
\begin{equation}
\bar{G}=\left[\bar{I}-\frac{1}{k^2}\nabla\nabla'\right]G\left(\bar{r},\bar{r}'\right)=\left[\bar{I}-\frac{1}{k^2}\nabla\nabla\right]\frac{e^{-jk|\bar{r}-\bar{r}'|}}{|\bar{r}-\bar{r}'|}
\end{equation}
子区域间的迭代公式为
\begin{equation}
\bar{T}_e\left(\bar{r},\bar{J}^{\left(k\right)}\right)=
-\bar{K}_e\left(\bar{r},\bar{J}^{\left(k-1\right)}\right)
+\left[\bar{E}^{inc}\left(\bar{r}\right)\right]_{tan},\bar{r}\in\Omega_i'
\end{equation}

对于MFIE有积分公式
\begin{equation}
\bar{K}_m\left(\bar{r},\bar{J}\right)=2\hat{n}\times\nabla\times\int_{\Omega_i}ds'G\left(\bar{r},\bar{r}'\right)\cdot\bar{J}\left(\bar{r}'\right),\bar{r}'\in\bar{\Omega_i},\bar{r}\in\Omega_i'
\end{equation}
\begin{equation}
\bar{T}_m\left(\bar{r},\bar{J}\right)=\bar{J}\left(\bar{r}\right)-2\hat{n}\times\nabla\times\int_{\Omega_i'}ds'G\left(\bar{r},\bar{r}'\right)\cdot\bar{J}\left(\bar{r}'\right),\bar{r}'\in\Omega_i',\bar{r}\in\Omega_i'
\end{equation}
同样，子区域间的迭代公式为
\begin{equation}
\bar{T}_m\left(\bar{r},\bar{J}^{\left(k\right)}\right)=
\bar{K}_m\left(\bar{r},\bar{J}^{\left(k-1\right)}\right)
+2\hat{n}\times\bar{H}^{inc}\left(\bar{r}\right),\bar{r}\in\Omega_i'
\end{equation}

CFIE的积分公式为
\begin{equation}
\bar{K}_c\left(\bar{r},\bar{J}\right)=\alpha\frac{\eta}{2}\bar{K}_m\left(\bar{r},\bar{J}\right)-\left(1-\alpha\right)\bar{K}_e\left(\bar{r},\bar{J}\right),\bar{J}\in\bar{\Omega_i},\bar{r}\in\Omega_i'
\end{equation}
\begin{equation}
\bar{T}_c\left(\bar{r},\bar{J}\right)=\alpha\frac{\eta}{2}\bar{T}_m\left(\bar{r},\bar{J}\right)+\left(1-\alpha\right)\bar{T}_e\left(\bar{r},\bar{J}\right),\bar{J}\in\Omega_i',\bar{r}\in\Omega_i'
\end{equation}
子区域间的迭代公式为
\begin{equation}
\bar{T}_c\left(\bar{r},\bar{J}^{\left(k\right)}\right)=
\bar{K}_c\left(\bar{r},\bar{J}^{\left(k-1\right)}\right)
+\alpha\eta\hat{n}\times\bar{H}^{inc}\left(\bar{r}\right)
+\left(1-\alpha\right)\left[\bar{E}^{inc}\left(\bar{r}\right)\right]_{tan},\bar{r}\in\Omega_i'
\end{equation}

重叠型区域分解法的求解步骤与经典的DDM是一致的，分为子区域内部求解的内迭代和子区域间相互影响的外迭代，具体实现如下：
\begin{enumerate}
\item 设各个扩展子区域的初始电流为零，外迭代次数$k=0$
\item $k=k+1$
\item for $i=1,2,\dots,M$，计算补余区域$\bar{\Omega}_i$上的电流对扩展子区域$\Omega_i'$上电流的作用，内部迭代求解$\Omega_i'$上的电流，最后舍弃缓冲区$\Omega_{bi}$上的电流，更新$\Omega_i$的电流
\item 计算相对余量误差，若误差达到精度要求，停止迭代，否则转到步骤2
\end{enumerate}

从数学角度来看，上述区域分解的迭代算法可以表示为矩阵公式
\begin{equation}\label{DDMiteration}
\tilde{Z}_{ii}\tilde{I}_i^{(k)}=\tilde{V}_i-\sum_{j<i,c(j)\notin b(i)}\tilde{Z}_{ij}I_j^{(k)}-\sum_{j>i,c(j)\notin b(i)}\tilde{Z}_{ij}I_j^{(k-1)},i=1,2,\dots,M
\end{equation}
其中，$b(i)$表示第$i$个子区$\Omega_i$缓冲区电流系数编号的集合，$c(j)$表示第$j$个子区$\Omega_j$上的电流系数，$\tilde{Z}_{ii}$与$\tilde{Z}_{ij}$由$\Omega_i$和缓冲区$\Omega_{bi}$的基函数，$\Omega_j$的基函数相互作用生成，分别为
\begin{equation}
\tilde{Z}_{ii}=\left[
\begin{array}{cc}
Z_{ii} & Z_{ib(i)}\\
Z_{b(i)i} & Z_{b(i)b(i)}
\end{array}
\right],\quad
\tilde{Z}_{ij}=\left[
\begin{array}{c}
Z_{ij}\\
Z_{b(i)j}
\end{array}
\right]
\end{equation}
式中，$Z_{ii},Z_{ij},Z_{b(i)b(i)},Z_{b(i)i},Z_{b(i)j}$分别是$\Omega_i$内基函数自作用，$\Omega_j$基函数对$\Omega_i$的作用，$\Omega_{bi}$内基函数自作用，$\Omega_i$基函数对$\Omega_{bi}$的作用，$\Omega_j$基函数对$\Omega_{bi}$的作用。
$I_j^{(k)}$是$\Omega_j$上第$k$次迭代求解的电流，$\tilde{I}_i^{(k)}$是$\Omega_i'$上第$k$次迭代求解的电流，$\tilde{V}_i$是$\Omega_i'$上的入射场。
它们分别表示为
\begin{equation}
\tilde{I}_i^{(k)}=\left[
\begin{array}{c}
I_i^{(k)}\\
I_{b(i)}
\end{array}
\right],\quad
\tilde{V}_i=\left[
\begin{array}{c}
V_i\\
V_{b(i)}
\end{array}
\right]
\end{equation}
式\ref{DDMiteration}实际上是变形的Gauss-Siedal迭代公式。
本文定义步骤4中的相对余量误差为电流的相对余量误差
\begin{equation}
error(k)=\frac{||I^{(k)}-I^{(k-1)}||}{||I^{(k)}||}
\end{equation}

由于区域分解算法是对上一次求得的电流进行修正，为了节省计算时间，可在第一轮外迭代时将各子区域内部的电流误差限放宽至$10^{-1}$，以后的每轮外迭代都缩小该误差限为上一轮的$\frac{1}{10}$，直到$10^{-3}$为止。
通常来说三次左右的外迭代即可使电流达到满意的精度范围内，因此该算法有着极好的计算效率。
\section{重叠型区域分解的缓冲区定义}
矩量法的离散化基于三角网格的剖分，那么基于矩量法的区域分解就必然涉及到将三角网格的分区和缓冲区的定义。

对于三角网格的分区可采用手动和自动两种方法：手动分区由操作人员指定区域的坐标范围，它具有划分区域边界规整的优点，但复杂目标的分区坐标定义比较繁琐，普适性不强；自动分区则是借助了开源软件gmsh的功能，它可以将已剖分好的网格按频谱、曲率、点的权重等各种算法进行分区，避免了手动指定坐标的过程，实现了对任意目标的自动化分区，有着较强的鲁棒性。

完成网格的分区后，需要对扩展缓冲区的方法进行定义，本文采用一类种子控制点的思想。
如下图\ref{figure:buffer_region}所示，粗线表示各子区域的边界，加重的点与子区域边界上的点同属一个三角形，但却不属于该子区域，这类点被称为子区域的“种子点”。
我们定义包含种子点却又不属于子区域的所有三角面片对为缓冲区，如图中的阴影部分所示。
\begin{figure}[htbp]
\centering
\includegraphics{buffer_region.jpg}
\caption{缓冲区定义}
\label{figure:buffer_region}
\end{figure}

上述一层种子点的定义还可以推广到两层甚至多层，但层数越多每个子区域的缓冲区就越大，求解的未知量就越多。
因此，通常对变化较平缓的区域边界只使用一层种子点，而边界电流变化较快的地方可以采用两层以上的种子点。

通过上述缓冲区的定义能够保证整个区域内部电流奇异性得到有效地抑制，与传统的非重叠型区域分解法相比，每个子区域只增加了少量的未知数，却使得算法的收敛性有了极大的提高。
从理论上讲，缓冲区越大子区域的计算量就越大，但收敛性却更好。
因此在实际问题中通常要折中考虑时间复杂度和空间复杂度，以达到最佳的计算效率。
\section{数值算例}
下面通过几个算例来验证算法的可靠性，方便起见，所有的算例都由$f=0.3$GHz的平面波照射，入射方向为$\theta_i=0^\circ,\phi_i=0^\circ$，观察点为$\phi=0^\circ,\theta=0$～$180^\circ$。
\subsection{开放结构的DDM}
算例1. 考虑一个边长为$2\lambda\times 3\lambda$的无限薄金属平板，平板上共有1839个未知数。
将平板按长边均匀分割成三块区域，模型和计算结果如图\ref{figure:DDM_plate}所示。
该算例的缓冲区仅在迭代子区域的前向与后向方向，可以证明本章所介绍的重叠型区域分解法为广义的前向后向缓冲迭代法(FBBS)，特殊情况下可自动退化为FBBS。
\begin{figure}[htbp]
\centering
\includegraphics[scale=0.5]{./data/chapter4/plate.jpg}
\includegraphics[scale=0.53]{./data/chapter4/plate.pdf}
\caption{边长$2\lambda\times 3\lambda$的平板}
\label{figure:DDM_plate}
\end{figure}

算例2. 半径为$1\lambda$的无限薄金属圆盘，共有933个未知量。
用相互垂直的两条直径将圆盘分为四块区域，在这种分区方式下各子区域就不存在前向后向的概念，此时FBBS方法将失去作用，但使用本章中重新定义的缓冲区扩展方式却依旧能处理该类问题。
模型和计算结果如图\ref{figure:DDM_disk}所示。
\begin{figure}[htbp]
\centering
\includegraphics[scale=0.57]{./data/chapter4/disk.jpg}
\includegraphics[scale=0.57]{./data/chapter4/disk.pdf}
\caption{半径为$1\lambda$的圆盘}
\label{figure:DDM_disk}
\end{figure}

算例3. 考虑非平面结构，一个边长为$1\lambda$的无盖立方体金属盒，未知量个数为1482。
金属盒壁为无限薄，按立方体的面将盖盒子划为五个区域，模型和计算结果如图\ref{figure:DDM_box}所示。
\begin{figure}[htbp]
\centering
\includegraphics[scale=0.5]{./data/chapter4/box.jpg}
\includegraphics[scale=0.6]{./data/chapter4/box.pdf}
\caption{边长为$1\lambda$的金属盒}
\label{figure:DDM_box}
\end{figure}
\subsection{封闭结构的DDM}
上述例子证明了重叠型区域分解法应用于开放结构物体中的正确性，下面将验证其在封闭结构物体中是否依旧有效。

算例4. 考虑一个金字塔型结构的金属散射体，塔底为边长$1\lambda\times 1\lambda$的矩形，塔尖距离塔底的垂直高度为$1\lambda$，整个物体的未知数为978。
此处将金字塔分为塔底和塔身两个区域，计算结果如图\ref{figure:DDM_pyramid}。
\begin{figure}[htbp]
\centering
\includegraphics[scale=0.5]{./data/chapter4/pyramid.jpg}
\includegraphics[scale=0.6]{./data/chapter4/pyramid.pdf}
\caption{底边长为$1\lambda$的金字塔}
\label{figure:DDM_pyramid}
\end{figure}

算例5. 柱型金属散射体，底面为半径$1\lambda$的圆，柱高$0.25\lambda$，共有2826个未知数。
按照划分金属圆盘的方式将其切成四块区域，模型和计算结果如图\ref{figure:DDM_cylinder}所示。
\begin{figure}[htbp]
\centering
\includegraphics[scale=0.53]{./data/chapter4/cylinder.jpg}
\includegraphics[scale=0.6]{./data/chapter4/cylinder.pdf}
\caption{底面半径为$1\lambda$的金属柱}
\label{figure:DDM_cylinder}
\end{figure}

\subsection{结果分析}
从上述的算例中可以看到，重叠型区域分解法的计算结果与原矩量法的结果吻合得很好，验证了其正确性，下面将就两种算法的效率进行对比。
\begin{table}[htbp]
\caption{DDM与普通矩量法的内存效率对比(MB)}
\label{DDM_memory}
\centering
\begin{tabular}{llllll}
\hline
算法 & 平板 & 圆盘 & 盒子 & 金字塔 & 圆柱\\
\hline
MoM & 54 & 14 & 35 & 15 & 126\\
DDM & 23 & 6 & 12 & 12.5 & 23\\
内存减少量 & 57.4\% & 57.1\% & 65.7\% & 20\% & 81.7\%\\
\hline
\end{tabular}
\end{table}

从表\ref{DDM_memory}中可以看到，使用DDM方法可以使内存的使用量得到极大的降低，这使得在同样的计算机硬件资源条件下能求解更大规模的电磁问题。
\begin{table}[htbp]
\caption{DDM与普通矩量法的时间效率对比(s)}
\centering
\begin{tabular}{llllll}
\hline
算法 & 平板 & 圆盘 & 盒子 & 金字塔 & 圆柱\\
\hline
MoM & 32.146 & 8.907 & 23.057 & 10.434 & 74.841\\
DDM & 12.917 & 3.272 & 16.825 & 27.521 & 72.949\\
时间减少量 & 59.8\% & 63.3\% & 27.0\% & -163.8\% & 2.5\%\\ 
\hline
\end{tabular}
\end{table}

而与内存使用量不同的是，计算效率的提升除了与子区域的大小相关外，还与模型形状和子区域边界所处的位置有关。
正如前文中提到的，平板、圆盘等结构的区域边界电流变化很小，收敛速度就很快，计算效率的提高很明显；盒子、圆柱等结构的区域边界电流变化较快，收敛速度较慢，计算效率的提高就很有限。
甚至在金字塔的例子上出现计算时间反而变长的情况，这是由于金字塔底和塔身的夹角较尖锐，区域边界的电流有较大突变造成的。
因此在选取区域边界的时候要尽量避开这些位置，通常可使用基于曲率半径的分区算法，以保证算法的计算效率。
\section{本章小结}
本章介绍了一种基于迭代求解的分治混合算法:区域分解法(DDM)。
首先通过数学公式严格证明和推导了矩阵分裂迭代求解的可行性和算法复杂度，而后将其思想应用到基于积分方程的矩量法中。
鉴于矩量法特殊的电流连续性要求，提出使用缓冲区对子区域进行扩展，以抑制边界上的电流奇异性，并简要阐述了扩展缓冲区的规则。
最后通过几个数值算例证明了该算法的正确性和效率，并指出了影响效率最大化的根本因素。
\chapter{基于等效源法的快速算法框架}
上一章节所介绍的区域分解法有着计算效率高、精度可靠和易于并行化等诸多优点。
通过区域分解，将整个求解空间上的未知电流划分到若干个耦合的子区域中，通过各子区域间解的迭代而最终达到稳态。

然而正如上章末尾处所述，区域分解算法的计算效率不仅和子区域的大小相关，还和子区域迭代矩阵的条件数紧密挂钩，而矩阵的条件数是影响其收敛速度的决定性因素。
通常来说，基于表面积分方程的算法中精度与收敛性总是一对不可调和的矛盾，很难同时达到最优。
为了解决上述问题，本章中将介绍一种基于Huygens原理的等效源法(EPA)，并在此算法的基础上建立一套快速算法的计算框架，将基于电流修正的混合方法与基于迭代的混合方法有机的结合起来。
\section{等效源法的原理与基本思想}
假设一组(P个)材质均匀的物体$D_p\left(p=1,2,\dots,P\right)$置于各向同性的均匀背景$D_0$中，背景和散射体的介电常数和磁导率分别为$\epsilon_p$和$\mu_p\left(p=0,1,\dots,P\right)$。
计算时谐场$\bar{E}_0^{inc}$和$\bar{H}_0^{inc}$入射情况下该场景的电磁散射问题。
假设入射波的时谐因子为$e^{j\omega t}$，其中$\omega=2\pi f$，$f$为电磁波的频率。

用$S_p$表示$D_p$的表面，由$D_p$产生的散射场可通过下列表达式得出
\begin{equation}
\bar{E}_p^{sca}\left(\bar{r}\right)=\eta_0\bar{\mathcal{L}}_p\left(\bar{J}_p\right)\left(\bar{r}\right)
-\bar{\mathcal{K}}_p\left(\bar{M}_p\right)\left(\bar{r}\right)
\end{equation}
\begin{equation}
\bar{H}_p^{sca}\left(\bar{r}\right)=\bar{\mathcal{L}}_p\left(\bar{M}_p\right)/\eta_0
+\bar{\mathcal{K}}_p\left(\bar{J}_p\right)\left(\bar{r}\right)
\end{equation}
其矩阵形式为
\begin{equation}\label{EPA:scattering}
\left[\begin{array}{c}
\bar{E}_p^{sca}\\
\eta_0\bar{H}_p^{sca}
\end{array}\right]=
\left[\begin{array}{cc}
\eta_0\bar{\mathcal{L}}_p & -\bar{\mathcal{K}}_p\\
\eta_0\bar{\mathcal{K}}_p & \bar{\mathcal{L}}_p
\end{array}\right]
\left[\begin{array}{c}
\bar{J}_p\\
\bar{M}_p
\end{array}\right]
\end{equation}
其中$\eta_0=\sqrt{\mu_0/\epsilon_0}$为$D_0$中的波阻抗，
\begin{equation}
\bar{J}_p=\hat{n}_p\times\bar{H}_p
\end{equation}
\begin{equation}
\bar{M}_p=-\hat{n}_p\times\bar{E}_p
\end{equation}
分别为$S_p$上感应电流和感应磁流的强度，$\hat{n}_p$是$S_p$上指向$D_0$的单位法向矢量。
面积分算子$\bar{\mathcal{L}}_p$和$\bar{\mathcal{K}}_p$定义如下
\begin{equation}
\begin{split}
\bar{\mathcal{L}}_p\left(\bar{F}\right)\left(\bar{r}\right)=&-\frac{j}{4\pi k_0}\nabla\int_{S_p}G_0\left(\bar{r},\bar{r}'\right)\nabla_S'\cdot\bar{F}\left(\bar{r}'\right)dS'\\
&-\frac{jk_0}{4\pi}\int_{S_p}G_0\left(\bar{r},\bar{r}'\right)\bar{F}\left(\bar{r}'\right)dS'
\end{split}
\end{equation}
\begin{equation}
\bar{\mathcal{K}}_p\left(\bar{F}\right)\left(\bar{r}\right)=\frac{1}{4\pi}\nabla\times\int_{S_p}G_0\left(\bar{r},\bar{r}'\right)\bar{F}\left(\bar{r}'\right)dS'
\end{equation}
这里的$\bar{F}$可以是$\bar{M}_p$或$\bar{J}_p$，$G_0$为环境$D_0$中的格林函数，$k_0=\omega\sqrt{\epsilon_0\mu_0}$为$D_0$中的波数。
因此，求$D_0$中的散射场的问题就转化为求解$S_p\left(p=1,2,\dots,P\right)$上的未知电流$\bar{J}_p$和未知磁流$\bar{M}_p$的问题。

根据等效原理，处于封闭区域内部或外部的场可以通过区域边界上场的切向分量决定。
因此，EPA的基本思想是将$S_p$上的未知量转移到虚拟表面$ES_p$上，以$ES_p$上的虚拟电磁流去等效原$S_p$上的感应电磁流，使得$ES_p$外部区域的散射场与原问题一致。

如图\ref{figure:EPA}所示，$ES_p$为单独包围$S_p$的等效面$\left(p=1,2,\dots,P\right)$,且不与任何$ES_l\left(l\ne p\right)$相交。
\begin{figure}[htbp]
\centering
\includegraphics[scale=0.2]{EPA.eps}
\caption{等效源法的基本思路}
\label{figure:EPA}
\end{figure}
在每个虚拟等效面ES上存在方程
\begin{equation}\label{EPA:equation}
\left[\begin{array}{c}
\bar{J}_p^{sca}\\
\bar{M}_p^{sca}
\end{array}\right]=
S_{pp}\left[\begin{array}{c}
\bar{J}_p^{inc}\\
\bar{M}_p^{inc}
\end{array}\right]+
S_{pp}\left(\sum_{l=1,l\ne p}^PT_{pl}\left[\begin{array}{c}
\bar{J}_l^{sca}\\
\bar{M}_l^{sca}
\end{array}\right]\right)
\end{equation}
其中$S_{pp}$表示等效面$ES_p$的自散射参数，即把$ES_p$上的入射源转换为散射源的算子；$T_{pl}$表示$ES_l$对$ES_p$的耦合参数，即把$ES_l$上的散射源转换为$ES_p$上入射源的算子。
式\ref{EPA:equation}中$\bar{J}^{inc}$和$\bar{M}^{inc}$为已知的入射源，因此只要知道参数$S$和$T$，就能求出ES上的虚拟感应流$\bar{J}^{sca}$和$\bar{M}^{sca}$，再通过\ref{EPA:scattering}式的形式便可求出$D_0$中的散射场。

求解感应流$\bar{J}^{sca}$和$\bar{M}^{sca}$的过程与传统矩量法类似，也是求解一个矩阵形式的线性方程组
\begin{equation}
\label{EPA:matrixequation}
\left[\begin{array}{cccc}
I & -S_{11}T_{12} & \dots & -S_{11}T_{1P}\\
-S_{22}T_{21} & I & \quad & -S_{22}T_{2P}\\
\vdots & \quad & \ddots & \vdots\\
-S_{11}T_{P1} & \dots & \quad & I
\end{array}\right]
\left[\begin{array}{c}
Q_1^{sca}\\
Q_2^{sca}\\
\vdots\\
Q_P^{sca}
\end{array}\right]
=\left[\begin{array}{c}
S_{11}Q_1^{inc}\\
S_{22}Q_2^{inc}\\
\vdots\\
S_{PP}Q_P^{inc}
\end{array}\right]
\end{equation}
这里的$I$是单位阵，
\begin{equation}
Q_p^{sca}=\left[\begin{array}{c}
\bar{J}_p^{sca}\\
\bar{M}_p^{sca}
\end{array}\right]=
\left[\begin{array}{c}
\hat{n}_p\times\bar{H}^{sca}\\
-\hat{n}_p\times\bar{E}^{sca}
\end{array}\right]
\end{equation}
和
\begin{equation}
Q_p^{inc}=\left[\begin{array}{c}
\bar{J}_p^{inc}\\
\bar{M}_p^{inc}
\end{array}\right]=
\left[\begin{array}{c}
\hat{n}_p\times\bar{H}^{inc}\\
-\hat{n}_p\times\bar{E}^{inc}
\end{array}\right]
\end{equation}
为$ES_p\left(p=1,2,\dots,P\right)$上的入射和散射等效流强度。
\section{自散射参数和互耦合参数}
EPA中定义自散射参数$S$可用方程式表示为
\begin{equation}\label{EPA:selfaffect}
S_{pp}=C_p\left(B_{pp}\right)^{-1}A_p
\end{equation}
其中，$A_p$为内向转移矩阵，将作用于$ES_p$的激励转换成$S_p$上的激励；$B_{pp}$为$S_p$的阻抗特性矩阵，其逆矩阵$\left(B_{pp}\right)^{-1}$将激励转换成$S_p$上的感应流；$C_p$为外向转移算子，将$S_p$上感应流转换成$ES_p$上散射源。
因此，$S$参数模拟了电磁波从外到内激励——在物体上产生感应流——从内到外散射这一物理过程。
它类似于微波网络中的自反射系数，因为对于$ES_p$的外部区域来说，我们只需要知道$S_{pp}$的表达式，而不需要关心$ES_p$内部具体是什么形式。

互耦合参数$T$用矩阵形式表达为
\begin{equation}\label{EPA:transfer}
T_{pl}=\left[\begin{array}{cc}
\hat{n}_p\times\bar{\mathcal{K}}_{pl} & \hat{n}_p\times\bar{\mathcal{L}}_{pl}/\eta_0\\
-\eta_0\hat{n}_p\times\bar{\mathcal{L}}_{pl} & \hat{n}_p\times\bar{\mathcal{K}}_{pl}
\end{array}\right]
\end{equation}
其中$\hat{n}_p$为$ES_p$的外向单位法向矢量，$T_{pl}$体现$ES_l$对$ES_p$的耦合作用。

上述\ref{EPA:selfaffect}、\ref{EPA:transfer}两式中，参数$T$的形式始终保持不变，参数$S$与物体材质有关。
当$S_p$为PEC时，
\begin{equation}
A_p=\left[\begin{array}{c}
a_p\left(\bar{\mathcal{L}}_p\right)_{tan}+b_p\hat{n}_p\times\bar{\mathcal{K}}_p\\
-a_p\left(\bar{\mathcal{K}}_p\right)_{tan}/\eta_0+b\hat{n}_p\times\bar{\mathcal{L}}_p/\eta_0
\end{array}\right]^T
\end{equation}
其中$a_p$和$b_p$分别为混合场积分方程(CFIE)中电场积分方程(EFIE)的系数和磁场积分方程(MFIE)的系数，$B_{pp}$为矩量法相应CFIE的系数矩阵$Z$，外向转移算子为
\begin{equation}
C_p=\left[\hat{n}_p\times\bar{\mathcal{K}}_p\quad -\eta_0\hat{n}_p\times\bar{\mathcal{L}}_p\right]^T
\end{equation}
\section{算法效率分析}
与传统方法相比，EPA多了$S$参数和$T$参数提取的前处理步骤，在计算时间上会略有增加，尤其在提取$S$参数时需要做矩阵的求逆操作，该操作的计算复杂度达到$O(N^3)$量级。
但是当求解的是由相同单元或者少数几类单元构成的大型周期或非周期结构问题时，每个基本单元的尺寸都比较小，且相同单元的逆矩阵只需求解一次，因此时间上所带来的额外开销也并不大。
重要的是，这些额外的开销可以带来矩阵条件数的极大改善，从而大幅减少问题的计算时间。
因为实际问题中常使用球形、柱形、立方体等经典结构做虚拟表面，这些结构的表面流变化都较平缓，可以进行粗剖分，且阻抗矩阵的条件数也比较好。
通常来说，限制大未知量问题求解速度的主要瓶颈在于矩阵的条件数太差，这导致其收敛速度慢，需要过多的迭代次数，这也正是为什么许多人致力于研究优秀的预处理器。
而通过使用EPA算法，可以做到在不降低剖分密度(剖分密度与计算精度紧密相关)的情况下改善矩阵的性态，有效地解决了可靠性与收敛性无法同时达到最优的问题，从而达到提高计算效率的目的。

同时，我们假设散射体表面$S_p$上有$N_p$个未知量$\left(p=1,2,\dots,P\right)$，其相应的等效面$ES_p$上有$M_p$个未知量。
由于$ES_p$处于原问题的近场区，均匀空间中的场值总是光滑连续的，所以在使用上文所述的几种虚拟表面结构时，$ES_p$上的网格可以不用严格遵循$0.1\lambda$的剖分准则，通常来说$M_p$远小于$N_p$。
基于表面积分方程的算法所需要的存储空间为$O(N^2)$量级($N$为未知量的个数)，因此通过使用EPA还可以达到减少未知量、从而降低内存需求的目的。
\section{EPA算法的应用拓展}
EPA算法除了上述改善矩阵条件数的作用外，还可以与前两章所阐述的算法相结合，延伸出更高效的解决方案。
\subsection{EPA与PO-MM结合}
第三章中提到，为了对矩量法区做PO电流修正，需要计算两个区域间的耦合矩阵，而PO-MM混合算法的最大瓶颈就在于这两个大型矩阵的乘法运算，这个过程耗费了过多的时间。
假设MM区有未知量$N$个，PO区有未知量$M$个，则这两个矩阵的规模分别为$N\times M$和$M\times N$，其相乘的计算量为$N^2M$，可见减少MM区的未知量数目可对算法进行平方加速。
而EPA算法恰好具有这个功能，因此可在每个矩量法区用等效面代替原物体，将PO区的效应耦合到ES上，并在ES上迭代求解该系统方程。
另外由于$N$的减小，等效面上$N\times N$的阻抗矩阵、两个区域间$N\times M$和$M\times N$的耦合矩阵的规模也都随之减小，因而同时降低了算法对内存的需求量。
\subsection{EPA与DDM结合}
原始的EPA是经由求解式\ref{EPA:matrixequation}的大型矩阵而得到空间$D_0$中的场。
然而通过观察式\ref{EPA:equation}的形式，发现其具有可区域分解的特性，因此结合上一章中的算法，可以进一步对EPA算法的效率进行优化：
只要提取出所有的$S_{pp}$和$T_{pl}(p=1,2,\dots,P,l\ne p)$参数后，便可使用DDM的思想迭代求解问题，具体步骤如图\ref{fig:EPA_flow_chart}所示。
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=1.0]{EPA_flow.eps}
  \caption{EPA算法流程图}
  \label{fig:EPA_flow_chart}
\end{figure}
通常来说，ES上的散射源经过3次外迭代更新即可达到很高的精度。
如上一节中所论述的，参数提取过程是EPA方法中消耗时间和存储量最大的环节，但不同的参数$S_{pp}$和$T_{pl}(p=1,2,\dots,P,l\ne p)$之间是相互独立，可通过并行处理技术进行加速。
另外，求解式\ref{EPA:equation}的步骤同样可以并行化，但由于各线程运算效率不同，收敛速度可能会有所下降，通常要多更新1～2次ES上的散射源才能达到串行处理时同样的精度。
\subsection{基于EPA的快速算法框架}
通过上面的论述，我们可以看到所有基于矩量法及其混合的快速算法都可以通过使用EPA进行优化，因此可以用EPA算法作为整个快速算法框架的基础，这种框架有如下的优点：
\begin{itemize}
\item 任何使用矩量法的区域都用一个虚拟等效面包裹起来，而后所有的不论是DDM还是PO-MM都在虚拟面之间或虚拟面和PO区之间进行。
\item 可以统一化外部求解的整个过程，而不用针对矩量法区到底是PEC还是介质而使用不同的迭代方程。
\item 推广到更远来说，虚拟表面内部甚至可以不用矩量法求解，例如计算非均匀介质的时候可以用上FDTD或FEM。
\end{itemize}
通过使用EPA，所有的数值算法都可以收敛到基于并行区域分解的矩量法，大大简化了通用型仿真软件的设计流程。
\section{数值算例}
\subsection{正确性验证}
算例1. 首先先验证EPA算法的正确性，一个边长为$1\lambda$的PEC正立方体处于自由空间中，入射波为频率0.3GHz的平面波，入射角为$\theta_i=0^\circ,\phi_i=0^\circ$。
用一个半径为$1\lambda$的球壳型等效面包裹该立方体，球心与立方体中心重合。
经过剖分后，立方体表面有1809条公共边，球壳上有480条公共边(图\ref{fig:EPA_cube})，从这里就可以看到，使用EPA算法至少可以降未知数降到一半以下。
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.4]{./data/chapter5/cube.jpg}
  \includegraphics[scale=0.9]{./data/chapter5/sphere_cube.jpg}
  \caption{边长为$1\lambda$的正立方体及球壳型等效面}
  \label{fig:EPA_cube}
\end{figure}

分别计算HH和VV极化下$\theta=0$～$180^\circ$，$\phi=0^\circ$的雷达散射截面(RCS)，并与矩量法进行对比，结果吻合良好，证明了EPA的正确性(图\ref{fig:EPA_verification})。
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.57]{./data/chapter5/cube_HH.pdf}
  \includegraphics[scale=0.57]{./data/chapter5/cube_VV.pdf}
  \caption{正立方体的EPA算法验证}
  \label{fig:EPA_verification}
\end{figure}
\subsection{周期结构问题中的应用}
算例2. 周期结构问题的计算是EPA算法应用最广泛的地方，用上述相同的立方体单元分别组成$1\times 2$，$2\times 2$，$2\times 2\times 2$的阵列，排列方式如图\ref{fig:EPA_2cube}，\ref{fig:EPA_4cube}，\ref{fig:EPA_8cube}所示，每个立方体都包裹于相同的球壳中。
计算VV极化下$\theta=0$～$180^\circ$，$\phi=0^\circ$的雷达散射截面，并与矩量法结果进行对比，可以看到，在计算周期结构目标时EPA算法在精度上有着很好的表现。
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.3]{./data/chapter5/2cube.jpg}
  \includegraphics[scale=0.57]{./data/chapter5/2cube_VV.pdf}
  \newline
  \caption{$1\times 2$的立方体阵列}
  \label{fig:EPA_2cube}
\end{figure}
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.35]{./data/chapter5/4cube.jpg}
  \includegraphics[scale=0.57]{./data/chapter5/4cube_VV.pdf}
  \newline
  \caption{$2\times 2$的立方体阵列}
  \label{fig:EPA_4cube}
\end{figure}
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.35]{./data/chapter5/8cube.jpg}
  \includegraphics[scale=0.57]{./data/chapter5/8cube_VV.pdf}
  \newline
  \caption{$2\times 2\times 2$的立方体阵列}
  \label{fig:EPA_8cube}
\end{figure}

对比算例2中三个问题的效率，由于阵列具有周期结构，计算时使用OpenMP的并行处理技术，并与成熟商业仿真软件FEKO的矩量法和快速多极子方法(均开启并行计算)对比内存和计算时间的消耗，结果如图\ref{fig:EPA_efficiency}。

可以看到，EPA算法无论在内存还是时间效率上都远远高于普通矩量法，与并行快速多极子相当。
需要注意的是内存效率中对比的是最高内存使用量，事实上EPA算法仅有在计算自散射参数时需要大量内存，这个过程仅占算法的一小部分且可以提前进行预处理，因此实际上的内存效率和快速多极子相当。
另外从图上还可以看出EPA的算法效率并不随周期数的增加而大幅降低，符合该算法适合计算周期结构问题的结论。
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.57]{./data/chapter5/memory.pdf}
  \includegraphics[scale=0.57]{./data/chapter5/time.pdf}
  \caption{EPA的算法效率}
  \label{fig:EPA_efficiency}
\end{figure}
\newline

算例3. EPA算法中一个虚拟表面所包含的目标单元数可以有多个，这种情况通常出现在单元过于密集，且每个单元的未知量数较少的情况下，如在Metamaterials中常见的SRR阵列结构（如图\ref{fig:SRR}）。
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.55]{./data/chapter5/SRR_array.pdf}
  \includegraphics[scale=0.35]{./data/chapter5/SRR_EPA_cube2.pdf}
  \caption{SRR阵列及其EPA单元}
  \label{fig:SRR}
\end{figure}
每个SRR单元上的未知数有1014个，通常来说3000未知数以下的问题使用LU求逆的效率还是很高的，因此我们在每个EPA虚拟面中包含三个SRR单元。
原问题共有$6\times 6\times 3$个SRR单元，73008个未知量，计算入射波垂直阵列面入射时的VV极化RCS，结果如图\ref{fig:SRR_results}。
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.7]{./data/chapter5/SRR_6_6_3_EPA.pdf}
  \caption{SRR阵列的计算结果}
  \label{fig:SRR_results}
\end{figure}

可以看到，EPA算法的结果在$\phi=40$～$140^\circ$范围内与精确解吻合的较好，而在这个范围外的精度就差强人意了。
这是由于在这个例子中我们使用了长方体型的虚拟表面，其棱边上电流变化较球形虚拟面来说要大，因此在应用边界条件求等效流时会产生数值上的误差。
通常来说，在问题允许的情况下尽量使用球壳型ES，这样得到的结果精度更高。
在本例中，由于SRR单元之间靠得太近，只能退而求其次，使用长方体型的虚拟表面。
\subsection{应用于PO-MM混合法}
算例4. 再验证以EPA算法作为基础框架对PO-MM算法进行加速的正确性及效率。
同样的立方体下方放置一边长$5\lambda\times 5\lambda$的金属平板，立方体处用EPA虚拟表面包裹，平板区用PO方法求解。
同第二章中原始的PO-MM混合方法对比如下图\ref{fig:EPA_PO}，可以看到其不但结果吻合的很好，计算效率上更有了质的飞跃(表\ref{tab:EPA_PO})。
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.37]{./data/chapter5/EPA_PO.jpg}
  \includegraphics[scale=0.57]{./data/chapter5/PO_VV.pdf}
  \caption{EPA应用于PO-MM混合法}
  \label{fig:EPA_PO}
\end{figure}
\begin{table}[htbp]
  \caption{PO-MM混合在使用EPA前后的效率对比}
  \centering
  \begin{tabular}{lll}
    \hline
    计算方法 & 时间消耗(s) & 内存消耗(MByte)\\
    \hline
    使用EPA前 & 266.016 & 880\\
    使用EPA后 & 41.085 & 350\\
    \hline
  \end{tabular}
  \label{tab:EPA_PO}
\end{table}

算例5. 用EPA算法为第三章中f16战斗机的例子加速。
用边长为1.2m的立方体虚拟表面包裹飞机上的天线，立方体上网格的最长边为$0.3\lambda$，共有219条公共边，而原天线上有3072个未知数。
对比各自使用PO-MM、EPA加速的PO-MM和MLFMA的RCS和计算效率，模型（图\ref{fig:EPA_f16}）和结果如下（图\ref{fig:EPA_f16_results}、表\ref{tab:EPA_f16}）。
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.7]{./data/chapter5/f16_EPA.pdf}
  \caption{f16战斗机中天线用EPA包裹}
  \label{fig:EPA_f16}
\end{figure}
\begin{figure}[htbp]
  \centering
  \includegraphics[scale=0.57]{./data/chapter5/f16_EPA_PO_MM_VV.pdf}
  \includegraphics[scale=0.57]{./data/chapter5/f16_EPA_MLFMA_VV.pdf}
  \caption{EPA应用于求解f16战斗机的例子}
  \label{fig:EPA_f16_results}
\end{figure}
\begin{table}[htbp]
  \caption{f16战斗机应用EPA前后的效率对比}
  \centering
  \begin{tabular}{lll}
    \hline
    计算方法 & 时间消耗(s) & 内存消耗(MByte)\\
    \hline
    使用EPA前 & 983.407 & 6000\\
    使用EPA后 & 92.52 & 900\\
    \hline
  \end{tabular}
  \label{tab:EPA_f16}
\end{table}

从图中可以看出，EPA算法可以以降低极小的精度为代价，大幅缩减PO-MM算法的计算时间和内存消耗，其计算结果与精确解的误差也在可接受的范围内。
\section{本章小结}
本章中介绍了一种基于Huygens等效面的等效原理法(EPA)，详细阐述了此种算法的思想和数学原理，包括等效流的表达形式，自作用算子和互作用算子的提取，着重说明了EPA在周期结构问题中的应用前景。
紧接着通过若干数值算例验证了其精度及效率上的可靠性，并在此算法的基础上将前两章的内容融合起来，建立了一套简洁、通用的快速算法框架。
通过使用这种框架，可以把各种不同的算法最后统一到基于矩量法的快速并行区域分解迭代中，简化了程序的设计，具有很高的实用价值。
\chapter{结束语}
\begin{enumerate}
\item 本文内容总结：

本文针对现如今计算电磁学发展中遇到的电大尺寸问题进行了研究，阐述了传统解决方案存在的不足，详细介绍了若干弥补这些不足之处的混合算法。
本文从三个方面展开研究：
\begin{itemize}
\item 介绍了基于高频物理光学近似的PO-MM混合算法，将电大尺寸物体的影响通过电流修正法耦合到矩量法区，大大减少了原问题的计算量。
\item 从计算机科学中的分治思想出发，研究了基于矩量法表面积分方程的区域分解法(DDM)，并针对边界电流奇异性问题研究了基于种子点的缓冲区定义方式。
\item 实现了一种基于Huygens原理的等效源法(EPA)，解决了表面积分方程中精度与收敛性不可同时达到最优的矛盾，并建立了一套基于EPA的快速算法框架，便于实现算法的软件化。
\end{itemize}
\item 有待进一步研究的工作：
  \begin{itemize}
  \item 弹跳射线法(SBR)是一种综合了基于电流和基于射线的高频方法，较PO来说有更高的精度，因此在做高低频混合时可用SBR替代掉PO，以得到更为精准的解。
  \item 基于时域积分方程的矩量法可一次性求得多个频点上的解，且具有天生的区域分解特性，其源的激励方式与真实情况也更为接近，因此值得做进一步的研究。
  \end{itemize}
\end{enumerate}
\end{Main} % 结束正文

\begin{Thanks}
在东南大学读研的三年，是我人生中至关重要的时光。在论文完成之际，首先我要衷心的感谢我的导师崔铁军教授。
崔老师严谨的治学态度，渊博的学识，勇于进取、善于开拓的精神给我留下了及其深刻的印象，时刻激励着我在学术上勇攀高峰。
特别是崔老师为人正直、平易近人、严于律己宽以待人的高尚品德，更是我今后人生道路上学习的楷模。

其次，我要感谢俞文明老师。俞老师在科研中一丝不苟、勤勤恳恳的态度为我树立了学习的榜样，尤其他对工作的热情无时无刻不在感染着我。
在此论文研究中，俞老师给了很多建设性的建议，令我受益匪浅。从俞老师身上我学到了什么叫精益求精，无论做任何事都要力求做到最好。

然后我要感谢实验室的老师和同学们，在工作和生活中，他们也给我很大的帮助。
感谢赵博师兄带我打羽毛球，感谢张俊师兄和陈刚师兄给我的指导和帮助，感谢贾晓倩师姐和何厅厅师兄对传承西电精神所做出的杰出贡献，感谢述润同学对知识的孜孜不倦，感谢周晓建设的绿色实验室环境，感谢苏璞同学时常带我出去装文艺，感谢游检卫师弟、陈林辉师弟、史珺师妹、王爱丽师妹和罗莹师妹给实验室带来的欢声笑语。没有你们，我这三年的研究生生活不会如此精彩，感谢你们！

我还要感谢我的室友孙引进、甘杰和张克迪，感谢他们在生活上对我的帮助，感谢他们提供的水果、夜宵、音乐还有卧谈会，感谢你们！

最后我要感谢我的妈妈，感谢她的养育之恩和一直以来对我的支持和信任，无论遇到什么困难，家庭永远是我最终的依靠！

感谢所有关心、支持和帮助过我的人，感谢你们！
\end{Thanks}

\bibliography{seuthesis}

\begin{Appendix}
\chapter{Maxwell方程组及其边界条件}\label{MaxwellEquation}
在介电常数和磁导率分别为$\epsilon$和$\mu$的均匀介质中，电场强度和磁场强度满足如下频域Maxwell方程
\begin{equation}
  \nabla\times\bar{E}=-\bar{M}-j\omega\mu\bar{H}
\end{equation}
\begin{equation}
  \nabla\times\bar{H}=\bar{J}+j\omega\epsilon\bar{E}
\end{equation}
\begin{equation}
  \nabla\cdot\bar{D}=\sigma_e
\end{equation}
\begin{equation}
  \nabla\cdot\bar{B}=\sigma_m
\end{equation}
其中$\bar{D}=\epsilon\bar{E}$，$\bar{B}=\mu\bar{H}$，上式中已约去了时谐因子$e^{j\omega t}$。
磁流$\bar{M}$和磁荷密度$\sigma_m$是物理上不存在的量，此处仅用于辅助求解辐射和散射问题。

不同介质区域的交界面上满足广义的电磁场边界条件
\begin{equation}
  -\hat{n}\times\left(\bar{E}_2-\bar{E}_1\right)=\bar{M}_s
\end{equation}
\begin{equation}
  \hat{n}\times\left(\bar{H}_2-\bar{H}_1\right)=\bar{J}_s
\end{equation}
\begin{equation}
  \hat{n}\cdot\left(\bar{D}_2-\bar{D}_1\right)=\sigma_e
\end{equation}
\begin{equation}
  \hat{n}\cdot\left(\bar{B}_2-\bar{B}_1\right)=\sigma_m
\end{equation}
其中$\hat{n}$为边界面上由区域2指向区域1的法向矢量。对于介质(区域2)与理想电导体(PEC，区域1)的交界面，边界条件变为
\begin{equation}
  -\hat{n}\times\bar{E}_2=0
\end{equation}
\begin{equation}
  \hat{n}\times\bar{H}_2=\bar{J}_s
\end{equation}
\begin{equation}
  \hat{n}\cdot\bar{D}_2=\sigma_e
\end{equation}
\begin{equation}
  \hat{n}\cdot\bar{B}_2=0
\end{equation}
同理，对于介质与理想磁导体(PMC)的交界面有
\begin{equation}
  -\hat{n}\times\bar{E}_2=\bar{M}_s
\end{equation}
\begin{equation}
  \hat{n}\times\bar{H}_2=0
\end{equation}
\begin{equation}
  \hat{n}\cdot\bar{D}_2=0
\end{equation}
\begin{equation}
  \hat{n}\cdot\bar{B}_2=\sigma_m
\end{equation}
\chapter{RWG基函数}\label{RWG}
RWG电流基函数是1982年有Rao,Wilton,Glisson等人提出的定义在一对相邻平面三角面片$T_n^{+},T_n^{-}$上的基函数，其定义为：
\begin{equation}\label{RWG:definition}
\bar{f}_n\left(\bar{r}\right)=
\begin{cases}
\frac{l_n}{2A_n^{+}}\bar{\rho}_n^{+}, & \bar{r}\in T_n^{+}\\
\frac{l_n}{2A_n^{-}}\bar{\rho}_n^{-}, & \bar{r}\in T_n^{-}\\
0, & \textrm{otherwise}
\end{cases}
\end{equation}
其中$l_n$是三角面片对$T_n^{\pm}$公共边的边长，$A_n^{\pm}$分别是面元$T_n^{\pm}$的面积，$\bar{\rho}_n^{\pm}$的定义如图\ref{figure:RWG}所示。
\begin{figure}[htbp]
\centering
\includegraphics{RWG.eps}
\caption{RWG基函数}
\label{figure:RWG}
\end{figure}

由式\ref{RWG:definition}可见，除了公共边外，面片对上的电流没有垂直于其边界的分量，因此边界上不存在线电荷。
而在每一条公共边上有，
\begin{equation}
\bar{f}_n^{+}\cdot\hat{h}_n^{+}=\frac{l_n}{2A_n^{+}}\bar{\rho}_n^{+}\cdot\hat{h}_n^{+}=\frac{l_n h_n^{+}}{2A_n^{+}}=1
\end{equation}
\begin{equation}
\bar{f}_n^{-}\cdot\hat{h}_n^{-}=\frac{l_n}{2A_n^{-}}\bar{\rho}_n^{-}\cdot\hat{h}_n^{-}=\frac{l_n h_n^{-}}{2A_n^{+}}=1
\end{equation}
其中，$\hat{h}_n^{\pm}$是在$T_n^{\pm}$上垂直于公共边的单位矢量，$h_n^{\pm}$是$T_n^{\pm}$在公共边上的高。
可见RWG基函数在公共边上满足电流连续性。

对$\bar{f}_n$求面散度有，
\begin{equation}\label{RWG:diverse}
\nabla_s\cdot\bar{f}_n\left(\bar{r}\right)=
\begin{cases}
\frac{l_n}{A_n^{+}}, & \bar{r}\in T_n^{+}\\
-\frac{l_n}{A_n^{-}}, & \bar{r}\in T_n^{-}\\
0, & \textrm{otherwise}
\end{cases}
\end{equation}
电流散度对应的物理含义是电荷量，由式\ref{RWG:diverse}可见，在一对三角面元上总的电荷量为零。

综合上述特点，RWG基函数几乎可以模拟绝大多数三维物体上的表面电流，因此得到了广泛的应用。
\chapter{线天线的矩量法}\label{LineMoM}
线天线的矩量法参考了R. F. Harrington的《Field Computation by Moment Methods》，这里假设线天线是一维结构，即天线截面半径远小于天线长度$\left(a<<L\right)$，且电流方向与轴线方向相同，即$\bar{J}\left(\bar{r}\right)=\hat{l}I\left(\bar{r}\right)$，$\nabla=\hat{l}\frac{\partial}{\partial l}$。

电场积分公式为
\begin{equation}
\left(j\omega\bar{A}+\nabla\Phi\right)_{tan}=\bar{E}^{inc}_{tan}
\end{equation}
其中的位函数
\begin{equation}
\bar{A}\left(\bar{r}\right)=\frac{\mu}{4\pi}\int_{L'}I\left(\bar{r}'\right)\cdot G\left(\bar{r},\bar{r}'\right)dl'
\end{equation}
\begin{equation}
\Phi\left(\bar{r}\right)=\frac{j}{4\pi\epsilon\omega}\int_{L'}\frac{\partial I\left(\bar{r}'\right)}{\partial l'}\cdot G\left(\bar{r},\bar{r}'\right)dl'
\end{equation}

将天线分成N段，每段上的电流可以视为常熟，如图\ref{figure:line_antenna}所示。
第n段由起点$n^-$到终点$n^+$，$\Delta l_n$代表第n段的长度。
\begin{figure}[htbp]
\centering
\includegraphics{line_antenna.eps}
\caption{任意形状的线天线}
\label{figure:line_antenna}
\end{figure}

将线上电流用脉冲基函数$P_n$展开，差分代替微分算子，并进行Galerkin检验可将积分方程化为矩阵$\left[Z\right]\left[I\right]=\left[V\right]$形式，其中
\begin{equation}
Z_{mn}=j\omega\mu\Delta l_n\cdot\Delta l_m\Psi\left(n,m\right)+\frac{1}{j\omega\epsilon}\left[\Psi\left(n^+,m^+\right)-\Psi\left(n^-,m^+\right)-\Psi\left(n^+,m^-\right)+\Psi\left(n^-,m^-\right)\right]
\end{equation}
$\Psi$的表达式如下
\begin{equation}
\Psi\left(m,n\right)=\begin{cases}
\frac{1}{4\pi\Delta l_n}\int_{n^-}^{n^+}\frac{e^{-jkR_m}}{R_m}dl\approx\frac{e^{-jkR_{mn}}}{4\pi R_{mn}} & m\ne n\\
\frac{1}{2\pi\Delta l_n}\log\left(\frac{\Delta l_n}{a}\right)-\frac{jk}{4\pi} & m=n
\end{cases}
\end{equation}
其中$R_{mn}=|\bar{r}_m-\bar{r}_n|$。

方程的右边项为
\begin{equation}
V_m=\int_{L}\bar{E}^{inc}\left(\bar{r}\right)\cdot P_m\left(\bar{r}\right)\hat{l}dl
\end{equation}
通常$V_m$为第$m$段的外加电压，可取
\begin{equation}
V_m=\left[\begin{array}{c}
0 \\
\vdots \\
V_i \\
\vdots \\
0
\end{array}\right]
\end{equation}
\end{Appendix}

\newpage
\printindex % 索引

%\begin{Resume}
%作者简介
%\end{Resume}

\backcover % 封底
\end{document}
